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CHAPTER  I 


EXECUTIVE  SUMMARY 


This  work  is  focused  on  the  development  of  computational  models  for  the  simulation  of  blast 
waves,  heterogeneous  combustion  and  their  interaction  with  particle  dispersion.  In  particu¬ 
lar,  an  Eulerian-Lagrangian  (EE/EL)  hybrid  approach  is  formulated.  In  this  approach,  the 
concept  is  to  employ  the  Eulerian  framework  for  the  dense  part  of  the  particles.  The  Eulerian 
approach  (EE)  is  suitable  for  the  simulation  of  dense  regimes  with  a  particular  advantage 
of  efhciently  simulating  a  large  number  of  particles.  This  is  made  possible  by  assuming  a 
body  of  particles  instead  of  individually  tracking  each  particle.  However,  this  assumption  is 
weakened  and  becomes  inaccurate  in  a  dilute  regime.  At  this  point,  the  Lagrangian  formu¬ 
lation  becomes  superior.  Such  transitions  from  dense  to  dilute  is  a  common  encounter  in  the 
dispersion  of  particles  due  to  blast  waves.  Our  hybrid  implementation  takes  into  account 
the  regimes  of  the  particles  being  simulated  and  employs  the  appropriate  formulation,  which 
provides  the  advantages  of  both  worlds  while  maintaining  the  required  level  of  accuracy. 

To  achieve  this  goal,  the  implementation  of  the  approach  is  first  validated  in  three  classic 
test  cases  with  experimental  or  theoretical  validation.  The  first  is  the  Sedov  case  where  the 
pressure  decay  and  blast  wave  front  are  validated  based  on  analytical  solutions.  In  this  test 
case,  the  Adaptive  Mesh  Refinement  (AMR)  capabilities  are  also  validated.  The  second  is 
the  Boiko  case  where  experimental  data  for  the  dispersion  of  a  clond  of  particles  in  a  confined 
tube  are  available.  The  good  agreement  of  the  results  from  these  simulations  validate  the 
EE,  EL  and  hybrid  EE/EL  approach.  The  third  validation  case  is  the  Rouge  case  which 
employs  a  similar  setup  to  Boiko’s.  These  results  validate  the  EE/EL  approach  as  well  as 
the  implementation  of  the  Discrete  Element  Method  (DEM). 

After  validating  the  numerical  approach  and  its  implementation,  several  cases  with  rele¬ 
vance  to  blast  wave  applications  have  been  performed  and  analyzed.  A  first  study  simulates 
the  propagation  of  a  blast  wave  due  to  a  spherical  charge  of  TNT  in  a  confined  room  with 
a  window  or  vent.  The  simulation  captures  the  various  wave  interactions  as  well  as  the  wall 
reflections  to  form  pressure  focusing.  The  results  also  demonstrate  the  shape  imprinting  on 
the  wave  front  from  the  initial  geometrical  and  structural  features  of  the  room  and  charge. 
Another  study  observes  the  qualitative  differences  in  the  detonation  of  a  cylindrical  charge 
with  two  configurations  of  detonators.  The  first  configuration  consists  of  a  detonator  along 
the  axis  of  the  cylindrical  charge.  In  this  manner,  the  detonation  wave  travels  in  a  radially 
symmetric  manner.  This  allows  the  assumption  of  using  a  detonation  profile  at  the  moment 
the  charge  has  been  consnmed  completely  to  initialize  the  propagation  of  the  blast  wave. 
The  second  configuration  sets  the  detonator  at  the  center  of  the  charge  which  causes  two 
detonation  waves  propagating  axially  from  the  center  to  the  extremities.  In  this  case,  the 
condensed  phase  propagation  is  required  to  be  resolved,  which  renders  the  simulation  costly. 
The  results  show  a  substantial  influence  of  the  initial  stages  on  the  longer  term  blast  wave 
shape  and  propagation. 

A  more  complex  application  that  combines  the  features  of  the  previous  simulations  is 
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the  propagation  of  a  blast  wave  in  rooms  connected  with  hallways.  In  this  case,  C4  and 
Nitro-methane  charges  are  used.  The  simulation  captures  the  wave  interactions,  reflections, 
and  pressure  focusing,  and  carbon  soot  is  tracked.  Comparisons  of  pressure  confinement 
with  one  room  and  rooms  with  hallways  depict  the  expected  pressure  relief.  Additionally, 
the  effect  of  the  charge  shape,  after-burning,  and  spore  aerosol  modeling  has  been  studied. 

Building  on  these  studies,  a  two-room  configuration  is  implemented  to  study  the  disper¬ 
sion  of  particles  due  to  a  blast  wave  in  a  confined  room.  The  two-room  comes  into  play 
because  of  a  vent  in  the  first  room  that  would  disperse  the  soot  particles  and  pressure  wave 
into  the  larger  room.  A  cloud  of  spore  particles  is  initialized  in  a  corner  of  the  smaller  room 
to  investigate  their  survivability.  Results  show  hydrodynamic  instabilities  due  to  RMI  and 
RTI  particularly  at  the  locations  where  the  reflected  pressure  waves  encounter  detonation 
products.  When  the  blast  wave  contacts  the  spore  particles,  the  spore  particle  cloud  is  de¬ 
formed  along  the  contact  surface  of  the  shock.  The  spores  are  driven  towards  the  wall,  and 
pushed  onto  the  corner  of  the  inner  room  with  relatively  higher  gaseous  temperature. 

In  another  study,  wave  interactions  due  to  multi-blasts  are  investigated.  These  also 
include  pressure  wave  interactions  with  a  fixed  wall.  Finally,  a  study  employing  the  transition 
of  a  shock  to  detonation  and  then  to  a  blast  wave  is  performed  to  study  the  influence  of 
micro-structural  changes  in  the  energetic  material  on  the  long  term  blast  wave.  Three 
configurations  are  studied  where  the  charge  constitutes  of  either  pure  energetic  material  or 
contains  crystals  with  a  polymer  binder.  Additionally,  the  shock  location  is  varied.  The 
results  show  a  substantial  influence  of  minute  changes  in  the  initial  micro-structure  on  the 
resulting  blast  wave  and  turbulent  mixing. 
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CHAPTER  II 


INTRODUCTION/MOTIVATIONS 


Deployment  of  explosives  is  ubiquitous  in  several  engineering  industries  encompassing  nu¬ 
merous  different  applications  ranging  from  mining  to  modern  warfare,  from  fire  quenching  in 
oil  field  to  agent-defeat  etc.  Although  explosions  have  been  widely  studied  by  the  research 
community  for  well  over  a  century,  many  phenomena  still  remain  to  be  investigated  in  order 
to  properly  understand  and  characterize  the  how-field  in  the  post-detonation  regime.  The 
complexity  of  post-detonation  how  field  aided  by  vigorous  turbulent  mixing  can  be  distinc¬ 
tive  due  to  the  presence  of  many  scales  of  motions,  flow  jetting,  high  level  of  shearing  motion, 
plume  surface  interactions  etc.  All  these  imply  that  an  initial  quiescent  ambient  can  become 
highly  turbulent  in  the  post-detonation  regime  which  in  turn  can  affect  the  final  after  burn¬ 
ing  significantly.  In  addition  to  this,  when  detonation  is  being  initiated  in  a  confined  space, 
it  can  generate  higher  level  of  pressure  load.  When  the  incident  pressure  wave  impinges  on  a 
structure  that  is  not  parallel  to  the  direction  of  the  wave  travel,  it  is  reflected  and  reinforced, 
producing  what  is  known  as  reflected  pressure.  The  reflected  pressure  is  the  force  to  which 
the  structure  ultimately  responds. 

The  problem  of  explosions  in  the  presence  of  particles  waxes  the  complexity  by  an  order 
of  magnitude  as  the  flow  now  involves  strong  discontinuities  such  as  shocks  and  contact 
surfaces,  as  well  as  smooth  flow  regions  such  as  shear  layers.  Also,  efficient  measures  have  to 
be  taken  to  track  these  particles  accurately  without  penalizing  the  computational  efiiciency. 
In  general  these  particles  can  be  reactive  e.g.  soot  particles,  bio-agents  etc.  which  can  be 
ignited  by  blast  wave  and  eventually  combusted  via  its  entrainment  into  flames.  In  addition 
to  this,  more  than  one  type  of  particles  may  be  present  in  the  domain  e.g.  bio-agents  and 
Aluminum,  soot  and  Aluminum  etc.  which  might  react  with  each  other  as  well.  Overall, 
there  exists  a  plethora  of  challenges  when  it  comes  to  model  explosions  with  particles.  The 
dispersion  of  these  particles  is  highly  transient  involving  wide  gamut  of  length  and  time  scales. 
The  detonation  generated  plume  can  impact  the  dispersion,  mixing  and  burn-out.  Due  to 
this,  in-homogeneous  mixing  resulting  in  plume-induced  instability  (i.e.  RM  instability) 
can  aggravate  burn-out  or  augment  dispersion.  Confined  space  adds  to  the  flow  complexity 
by  introducing  corner  recirculation,  pressure  wave  reflections  etc.  Also,  in  confined  spaces, 
mixing  process  is  much  more  rapid  to  allow  combustion-induced  pressure  waves  to  coalesce 
and  create  sustained  pressure  loads  on  the  walls.  Characterization  of  the  particle  behavior 
in  the  post-detonation  regime  inside  a  confined  space  is  an  active  area  of  research  which 
requires  accounting  for  the  particle-turbulence-  shock/detonation-chemistry  interactions  in 
complex  three  dimensional  flows.  The  present  research  effort  aims  to  address  some  of  these 
issues. 

Experimental  measurements  can  be  highly  expensive  and  difficult  to  perform  as  the  am¬ 
bient  is  hostile.  Simulation  can  play  a  critical  role  in  this  scenario.  However,  simulating 
this  kind  of  complex  flow  mandates  to  use  higher  order  spatial  and  temporal  accuracy  that 
can  be  achieved  by  direct  numerical  simulation  (DNS)  or  large-eddy  simulation  (LES).  DNS 
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with  its  huge  computational  requirement  is  beyond  the  scope  of  current  and  perhaps,  the 
near-future  hardware  capability.  On  the  other  hand,  LES  of  these  complex  flows  requires 
advanced  sub-grid  models  as  well  as  a  robust  parallel  code.  Application  and  validation  of  the 
complex  turbulent  reacting  flows  have  been  demonstrated  previously  by  the  3D  LES  code 
developed  and  maintained  in  the  Computational  Combustion  Lab  (CCL)  at  Georgia  Insti¬ 
tute  of  Technology  (called  LESLIE3D,  hereafter)  under  DTRA  funded  work.  The  current 
effort  aims  to  evaluate  the  capability  of  the  code  by  validating  it  against  the  confined  space 
explosion  results.  The  next  step  is  to  enhance  the  capability  of  the  code  by  incorporating 
improved  sub-grid  models  for  calculating  particle  dispersion. 

The  goal  of  this  report  is  to  summarize  the  progress  made  so-far  in  order  to  build  a 
simulation  capability  to  perform  and  predict  the  complex  physics  of  detonations  and  blast 
waves  together  with  its  interaction  with  particle  dispersion  and  heterogeneous  combustion 
in  confined  spaces.  The  report  is  organized  as  follows:  Section  2.1.2  presents  the  technical 
objectives  and  statement  of  the  work.  Section  2.1.3  states  the  experimental  configuration 
which  is  followed  by  computational  setup  in  Section  2.1.4.  Section  2.1.5  deals  with  the 
numerical  setup.  Section  2.1.6  discusses  the  results.  The  report  is  wrapped  up  with  Section 
2.1.7  which  summarizes  future  plans  of  this  project. 
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CHAPTER  III 


OBJECTIVES 


The  primary  objective  of  this  effort  was  to  develop  and  deploy  a  state-of-the-art  and  next- 
generation  simulation  tool  for  WMD  studies  that  includes  advanced  subgrid  models  for 
turbulent  momentum  and  scalar  mixing,  and  combustion  in  complex  domains. 

The  proposed  study  will  leverage  the  availability  of  both  turbulent  mixing  models  and 
ANN  based  reaction  kinetics  and  enhance  them  with  new  subgrid  models  for  spores.  Both 
detonations  and  shear  turbulence,  are  captured  in  the  solver  using  a  hybrid  approach  that 
combines  a  robust  shock  capturing  model  with  a  less  dissipative  0(4)  central  scheme  in 
regions  without  strong  discontinuities.  A  Lagrangian  particle-tracking  model  can  track  liquid 
and/or  solid  particles  with  full  coupling  between  the  phases  [25,  4,  6].  This  solver  has  been 
successfully  employed  to  simulate  RM  instability,  gas  and  two- phase  detonations.  Effect  of 
collisions  between  parcels  and  within  the  parcel  is  included.  This  method  will  allow  modeling 
the  effect  of  collisions,  non-spherical  particles  and  local  clustering. 

Whereas  the  research  in  Georgia  Tech  is  basic  (unclassified)  research  investigating  fun¬ 
damental  mixing  and  combustion  processes  our  partner  CRAFT-Tech  has  developed  an  ad¬ 
vanced  CFD  code  capable  of  simulating  many  complex  (and  restricted)  problems  of  interest 
to  DTRA,  such  as  tests  at  Eglin  AFB  and  Ladeburg  Bunker  in  Germany.  The  GRAFT- 
Tech  code  is  well  established  with  a  long  history  of  application  to  many  classified  test  cases. 
This  combined  effort  will  add  sub  models  developed  in  Georgia  Tech  to  enhance  this  codes 
capability  for  these  applications. 

The  following  technical  objectives  apply  to  both  members  of  this  research  team. 

1.  Establish  a  new  coupled  Eulerian-Lagrangian  approach  to  track  dense  clouds  of  parti¬ 
cles  with  particular  emphasis  on  subgrid  models  for  spore  physics,  spore  distribution 
and  its  response  to  turbulent  mixing  in  the  post  detonation  flow  field. 

A  dense  Lagrangian  particle-tracking  model  is  being  validated  in  Georgia  Tech  cur¬ 
rently  under  exiting  DTRA  funding  [15,  16].  This  method  track  parcels  as  in  a  La¬ 
grangian  approach  but  also  includes  a  subgrid  stochastic  model  to  account  for  local 
clustering.  A  key  advantage  of  this  model  is  that  it  can  approximate  the  presence 
of  large  number  of  small  particles  in  a  subgrid  cloud  without  explicitly  tracking  each 
of  the  particles.  Thus,  it  has  the  ability  to  track  a  very  large  number  of  particles  in 
a  deterministic  but  cost-effective  manner.  This  approach  offers  an  approach  to  track 
particle  distribution  changes  locally  (due  to  afterburning,  collisions/merger,  breakup, 
etc.)  unlike  a  classical  volume  fraction  based  Eulerian  approach.  Therefore,  it  has  the 
potential  to  deal  with  spore  clouds  as  well  as  A1  particles  in  the  explosives  and  their  in¬ 
teractions.  The  current  approach  is  very  fundamental  (studies  in  isotropic  turbulence, 
spray  dispersion  etc)  but  will  have  to  be  generalized  to  deal  with  walls  etc.  Extension 
of  the  model  to  handle  subgrid  turbulence  (i.e.,  turbulence  at  the  small  scales),  spore 
sub-models  (including  kinetics) 
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2.  Validate  and  implement  computationally  elRcient  rate  kinetics  using  advanced  lookup 
approaches  based  on  ANNs  for  thermobaric  explosive  and  spore  mechanism. 

To  study  the  physics  of  thermobaric  explosives,  detailed  (or  at  least  multi-step)  kinetics 
or  finite  rate  models  are  needed  for  accurate  prediction  of  afterburning,  spore  thermal 
and  chemical  neutralization  mechanisms.  Due  to  a  wide  range  of  scales  computational 
cost  of  full-scale  problems  can  become  prohibitively  expensive,  and  on  top  of  this 
there  is  a  need  for  proper  treatment  under  turbulent  mixing  conditions.  The  subgrid 
model  developed  in  Georgia  Tech  has  combined  a  turbulent  mixing  model  with  finite 
rate  kinetics  and  also  developed  efficient  ANNs  to  paratemetrize  rate  mechanisms. 
This  strategy  will  be  extended  for  the  current  applications  and  integrated  into  the 
simulation  codes. 

3.  Integration  of  turbulent  mixing  and  combustion  models  into  full-scale  production  code 
for  DTRA  related  restricted  application  studies  This  objective  will  transition  all  the 
models  developed  during  this  research  into  the  CRAFT  Tech  production  code  and  then 
this  will  be  used  to  study  various  test  cases  identified  by  DTRA.  Further  refinements 
may  be  needed  and  will  be  part  of  this  overall  objective.  Another  part  of  this  effort 
will  be  the  release  of  some  versions  of  the  codes  for  DTRA  use  in  Federal  Laboratories, 
as  requested. 

The  goals  are  specifically  driven  by  DTRA’s  requirement  to  obtain  greater  confidence 
in  the  accuracy  of  computational  simulations  involving  advanced  high  energy  explosives  in 
closed  compartments/structures  and  agent  defeat  simulations  with  lower  energy  release  and 
dispersion  for  neutralizers.  The  overarching  technical  objective  is  to  assess  the  ability  in 
capturing  the  effects  of  turbulent  mixing,  disparate  time-scales,  multi-phase  physics  and 
scenario  input  uncertainty  on  AD  in  a  complex  geometry. 
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CHAPTER  IV 


MATHEMATICAL  FORMULATION  AND 
NUMERICAL  METHODOLOGY 


4.1  Modeling 

There  are  different  modeling  approaches  for  simulating  two-phase  flows,  which  depend  upon 
parameters  such  as  mass  and  volume  loading,  Stokes  number,  spray  regime  and  the  method 
used  for  simulating  the  carrier  phase  [37,  38,  39].  Some  of  the  typical  approaches  to  numer¬ 
ically  model  practical  multiphase  flow  systems  are  denoted  as  Euler ian-Lagrangian  (EL), 
Eulerian-Eulerian  (EE),  hybrid  and  statistical  methods  (others  may  exist  but  are  not  ad¬ 
dressed  here).  Table  4.1  summarizes  some  of  these  models  used  to  study  swirling  spray 
combustion  (non-swirling  spray  mixing  and  combustion  are  not  directly  addressed  here).  A 
brief  description  of  these  methods  with  their  advantages  and  limitations  is  provided  below. 

In  the  Eulerian-Eulerian  approach,  also  referred  to  as  the  two-fluid  method  [40,  41,  38, 
39],  both  carrier  and  disperse  phases  are  solved  using  an  Eulerian  framework.  The  two 
phases  are  considered  to  be  interpenetrating  and  requires  transport  equations  for  the  volume 
fraction,  velocity,  temperature  and  moments  of  size  distribution  of  the  dispersed  phase, 
which  is  obtained  after  homogenization.  The  mass,  momentum  and  energy  exchange  with 
the  carrier  phase  occurs  through  source/sink  terms.  Since  the  method  utilizes  a  common 
Eulerian  framework,  consistent  numerical  method  can  be  used  for  both  phases,  thus  leading 
to  easy  implementation  and  scalable  high-performance  parallel  computing.  However,  the 
method  requires  substantial  initial  modeling  effort  to  obtain  the  aforementioned  transport 
equations  for  the  dispersed  phase.  Additionally,  the  method  is  considered  to  be  expensive 
for  polydisperse  systems,  although  there  have  been  some  development  to  deal  with  multisize 
particle  sprays  in  a  computationally  efficient  manner  [42] .  Another  limitation  of  the  method 
is  associated  with  the  numerical  stability,  which  requires  that  the  concentration  gradient 
should  not  be  very  high  for  the  disperse  phase  within  the  flow  system. 

In  the  Eulerian-Lagrangian  method,  Lagrangian  tracking  of  the  dispersed  phase  is  per¬ 
formed  and  the  carrier  phase  is  simulated  using  the  conventional  Eulerian  framework  [43, 
44,  45,  46].  The  EL  method  is  the  most  common  approach  to  simulate  flow  systems  con¬ 
sidered  here,  due  to  its  robustness,  accuracy  and  ability  to  model  complex  phenomena  such 
as  poly-dispersity,  particle/wall  and  particle-particle  interactions,  crossing  trajectories,  and 
particle  break-up.  However,  the  method  is  computationally  expensive,  as  large  number  of 
particles  are  required  within  each  computational  cell  used  by  the  carrier  phase  to  allow  for 
a  smooth  Eulerian  reconstruction  of  the  feedback  force  to  the  carrier  phase. 

Another  major  issue  with  the  EL  method  is  related  to  the  accuracy  of  the  dispersed 
phase  statistics,  which  is  dependent  on  the  number  of  particles,  particularly  in  regions  of  a 
physically  observed  sparse  distribution  of  the  disperse  phase.  In  general,  the  convergence 
rate  of  the  dispersed  phase  statistics  by  the  EL  method  scales  as  where  N  is  the 
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Table  4.1:  Modeling  approaches  for  simulating  typical  swirling  spray  combustors. 


Method 

References 

Remarks 

Eulerian-Eulerian 

[4,  5,  6,  7,  8,  9,  10,  11,  12] 

Scalable  parallelization, 
consistent  numerical  method, 
subtantial  modeling  effort, 
expensive  for  polydisperse  sys¬ 
tems 

Eulerian- 

[13,  14,  15,  16,  17,  18,  19, 

20,  21] 

Robust  for  complex  systems, 

Lagrangian 

[22,  8,  10,  23,  24,  25,  26,  27 

,  28,  29] 

slower  statistical  convergence, 
inefficient  parallelization. 

Hybrid 

[30,  31] 

computationally  expensive 

Statistical 

[32,  33,  34,  35,  36] 

Consistent  numerical  method, 
closure  is  simpler, 
easy  to  parallelize, 
computationally  expensive 

number  of  particles  in  any  parcel  of  the  domain.  Regions  that  have  low  number  density  tend 
to  be  poorly  converged  because  N  is  small  in  comparison  with  regions  where  N  is  large. 

The  density-weighted  spatially  filtered  LES  equations  can  be  obtained  from  the  compress¬ 
ible  form  of  the  multi-species  Navier-Stokes  equations  by  employing  the  well  known  Favre 
filtering  approach.  The  Favre  filtered  (resolved)  quantity  corresponding  to  a  field  variable 
is  defined  as 

(f){x,t)  =  — - -  [  G{x,x' —  x)p{x' ,t)(f){x' ,t)dx' ,  (4.1) 

t)  Jo. 

where  Q,  is  the  computational  domain,  (S'  is  a  spatial  filter  function,  p  is  the  density  and 
(.)  denotes  the  conventional  spatial  filtering,  which  when  applied  to  the  density  field  p{x,t) 
leads  to  the  spatially  filtered  density  field  p{x,t)  given  by 

p{x,t)  —  /  G{x,x'  —  x)p{x' ,t)  dx' .  (4.2) 

Jn 

Note  that  the  two  filtering  approaches  are  related  through  pcj)  =  pcp. 

The  interphase  coupling  of  the  gas  phase,  the  Lagrangian  dispersed  phase  and  the  Eule- 
rian  dispersed  phase  can  be  tracked  in  terms  of  the  volume  fraction  through 

ag  -\-  ap  -\-  ad  =  1,  (4-3) 

where  ag,  cXp  and  Txg  denote  volume  fraction  of  the  gas  phase,  Lagrangian  dispersed  phase 
and  Eulerian  dispersed  phase,  respectively.  The  volume  fraction  of  the  Lagrangian  dispersed 
phase  cip  is  defined  as 
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(4.4) 


where  AV  is  the  volume  of  the  computational  cell,  Up  is  the  number  of  particles  in  each 
particle  group  (referred  to  as  a  parcel)  and  N  is  the  total  number  of  parcels  in  the  compu¬ 
tational  cell.  The  volume  fraction  of  the  gas  phase  ag  evolves  in  space  and  time  according 
to 


dag 

dt 


+  'W/,i 


dag 

dxi 


0, 


(4.5) 


where  is  the  interface  velocity  held.  From  here  onwards,  subscripts  'g',  ‘p’,  ‘d’  and  ‘7’  are 
used  to  indicate  the  gas  phase,  Lagrangian  particle,  Eulerian  dispersed  phase  and  interface 
quantity,  respectively.  As  mentioned  before,  the  governing  system  of  equations  for  EE  and 
EL  formulations  can  be  recovered  by  setting  =  0  and  =  0,  respectively,  in  the  general 
formulation  presented  below.  Note  that,  typically  for  EL  formulation  for  systems  having  low 
mass  or  volume  loading,  =  1  is  used. 


4.1.1  Eulerian  Gas  Phase 


Applying  a  box  hlter  (appropriate  for  a  hnite-volume  based  implementation)  to  the  system 
of  transport  equations  yields  the  density-weighted  hltered  LES  equations  for  the  gas  phase, 
which  comprises  of  transport  equations  for  mass,  momentum,  energy  and  species  expressed 
as 


dag  Pg  dag  pgltg^l  — 

+  dXi  “ 


dt 


dxj 

_  dag  —  dag  • 


dOLg  PqEg  d 


^9 

dt 


+ 


dXj 


Q!,q 


Pg^gJ^g  +  ^gjPg  +  Qgj  ^g^i^gji  + 


~  -  7-  -T- 


dag  PgYg  If  d 

dxi 


^9  g 

dt 


[Pg(yg.ku,,  +  Yg^,Vg,,,)  +  yj*  + 


agiOg  j^  -\-  for  k  1,2,...,  Ng, 


(4.6a) 


(4.6b) 


(4.6c) 


(4.6d) 


where  subscript  ‘71’  denote  the  source  term  contribution  from  Lagrangian  and  Eulerian 
dispersed  phases  and  Ng  is  the  number  of  species.  In  Eq.  (4.6)  Ui  is  the  velocity  vector, 
p  is  the  pressure,  E  is  the  total  energy,  is  the  mass  fraction  of  the  species,  is 
the  viscous  stress  tensor,  g,  is  the  heat-flux  vector  and  5ij  is  the  Kronecker  delta.  The 
subgrid  scale  terms  are  denoted  by  superscript  'sgs^  in  Eq.  (4.6)  and  they  require  closure 
approximation.  These  terms  include  subgrid  stress  tensor  r*?*,  subgrid  enthalpy  flux 
subgrid  viscous  work  o’*®*,  subgrid  convective  species  flux  Y)*®*  and  subgrid  diffusive  species 
flux  0®^®.  The  closure  models  for  these  terms  are  discussed  in  sec.  4.3.  The  source  terms  pp, 

F D,ii  Wd  and  SD,ki  which  appear  in  the  transport  equations  for  mass,  momentum,  energy 
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and  species,  respectively,  represent  the  interphase  exchange  terms.  They  are  obtained  in 
terms  of  the  source  terms  for  Lagrangian  and  Eulerian  dispersed  phases  through 


T 


D.m 


—  ^ p,m  +  -T , 


d.m  5 


(4.7) 


where  the  vector  =  [p,  ,  Q,  W,  represents  source  term  in  different  transport  equa¬ 
tions.  Here,  F  =  [Fi,  F2, . . . ,  and  S  =  [Si,  S2,  ■  ■ . ,  with  ‘dim’  denoting  the 

dimension  of  the  considered  flow  system  (2D/3D).  Note  that  for  mass  transfer  of  single 

species,  we  obtain  Su  =  Pd-  The  contribution  to  source  terms  from  the  Lagrangian  and 
Eulerian  dispersed  phases  are  provided  in  Eqs.  (4.15)  and  (4.14),  respectively. 

By  assuming  a  Newtonian  fluid  with  the  Stokes’  hypothesis  and  the  Fourier’s  law  of 
thermal  conduction,  the  Altered  viscous  stress  tensor  Tg^ij  and  the  heat  flux  vector  can 
be  approximated  as 


=  -^9 


dUg^i  , 

dxi 


+ 


9,9 


dxi 


2  dUg]^ 

3^^  dxk 


dXj 


Ns 


Yg,khg,kVgj^k  + 


(4.8a) 

SOS 

(4.8b) 

where  Pg  =  Pg{Tg)  is  the  dynamic  viscosity  of  the  gas.  Kg  is  the  thermal  conductivity  of  the 
gas,  hg^k  is  the  resolved  specific  enthalpy  of  the  species,  is  the  subgrid  heat  flux 

vector  (described  in  sec.  4.3)  and  i®  resolved  species  diffusion  velocities  that  can  be 
modeled  through  a  Fickian  diffusion  approximation  through 


^9,j,k  — 


Dg,k  dYg^k 

Y,.k  ' 


(4.9) 


Here,  Dg^k  denotes  diffusion  coefficient  of  the  species  and  it  can  be  obtained  from  a 
constant  Lewis  number  (Le)  assumption. 


4.1.2  Lagrangian  Particle  Phase 


The  Lagrangian  equations  for  the  motion  of  a  single  particle  within  dense  spray  combustors 
are  obtained  with  the  assumption  that  the  particle  density  is  much  greater  than  the  carrier 
(gas)  phase  density  {Pp/ pg  ~  10^)  and  particle  diameter  is  smaller  than  the  Kolmogorov 
length  scale.  These  equations  are  given  by 


^^p,i 

NT  ^ 

dnip  d  f  4:  o 


m. 


dUp  I  TT  9 y  {■ — '  //  I  /' — '  //  \  ^  a 

~  ^DPg\^g,i  4“  ^g,i  I  y^g,i  T  ^g^i  ^P,i)  3^^^  Qx  ^^p-^c,i-i 


dt 

dTj 


Tn„Cri—r^  =  2TrrpKgNu  (Tg  —  Tp)  —  ifipLy, 


(4.10a) 

(4.10b) 

(4.10c) 

(4.10d) 
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where  Xj,  m,  r,  T  denote  position  vector,  mass,  radius  and  temperature,  respectively.  Some 
of  the  other  terms  in  Eq.  (4.10)  include  the  drag  coefficient  Cd,  acceleration  due  to  inter¬ 
particle  interactions  the  heat  capacity  Cp,  the  thermal  conductivity  of  the  gas  phase 
Hg,  the  Nusselt  number  Nu  and  the  latent  heat  of  vaporization  L„.  The  Nusselt  number 
and  the  drag  coefficient  are  typically  expressed  as  empirical  functions  of  Reynolds  number 
(Re),  Prandtl  number  (Rr),  Mach  number  (M)  and  volume  fraction  ag  [47,  48,  49].  The 
closures  for  interface  quantities  can  be  found  elsewhere  [47,  48,  50].  The  acceleration  due 
inter-particle  interactions,  is  computed  as 


1  dr 

^pPp  dXi 


where  r  is  the  inter-granular  stress  given  by 

P  a  ^ 

r  =  - , 

^CS  ^p 


(4.11) 


(4.12) 


with  acs  being  the  close  packing  volume  fraction.  Also,  Pg  and  /5  are  empirical  constants, 
which  are  closed  based  on  the  nature  of  the  flow  being  considered  [51,  48,  50,  47]. 

In  Eq.  (4.10)(c),  the  sum  (ug^i  +  u”,^)  represents  instantaneous  (ug^i)  gas-phase  velocity 
components,  consisting  of  both  the  LES  resolved  velocity  Ug^i  and  the  unresolved  fluctuating 
velocity  uT,  which  can  be  reconstructed  by  employing  a  stochastic  model  [52,  23,  34].  For 
example,  the  unresolved  velocity  held  can  be  obtained  from  the  subgrid-scale  turbulent 
kinetic  energy  at  intervals  coincident  with  the  local  characteristic  eddy  lifetime  [52]  or  it  can 
be  based  on  a  stochastic  Markovian  model  [53] .  In  several  other  numerical  implementations  in 
the  past,  this  term  has  been  ignored,  however,  as  mentioned  in  [23],  in  poorly  resolved  regions 
of  the  flow,  where  the  subgrid-scale  turbulent  kinetic  energy  is  more  than  30%  of  the  resolved 
turbulent  kinetic  energy,  the  effect  of  unresolved  velocity  fluctuations  on  the  particle  motion 
becomes  important.  In  a  similar  way,  the  gas  phase  temperature  at  the  particle  location  in 
Eq.  (4.10)(d),  ideally  should  include  the  unresolved  temperature  fluctuation,  which  in  turn 
can  be  modeled  through  stochastic  means  [34] ,  however,  in  the  present  formulation,  we  have 
ignored  such  contribution. 


4.1.3  Eulerian  Particle  Phase 

In  the  limit  of  a  dense  mass  or  volume  loading,  the  dispersed  phase  can  be  modeled  as 
a  continuum  fluid.  In  such  cases,  the  governing  equations  for  the  dispersed  phase  can  be 
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expressed  in  the  Eulerian  framework  and  are  given  by 


doid  ^  dcx^  p^u^  i  _ 

Pd  1 


dt 


dxi 


doid  PijUd  i  ^  \ ~  ~  I  —  X  I  ^9^  —  M 

- Ft - ^  dFj  yPd^d,iUd,j  +  PdOij  +  -  Td,ij)\ 

_  Ocxd  „  _  doid  • 

dad  PdEd  ,  f  -  ~ 

dt  dxj 


(pd‘^d,jEd  +  Ud,jPd  +  qd,j  ~  UddTd,ji  +  +  <^2 


sgs 

j 


=  PlUlJ 


dad 

dxj 


Ul,iT  idj 


dad 

dxj 


Qd  -  ^d, 


with  the  source  terms  in  the  above  equations  given  by 

d  /4 


Pd  =  Ndirid  =  Nd—  i^Trrlpa  )  , 


Fdd  —  Nd 


TT 


rridUdd  +  ■^rjCDpg\ud,i  -  Ug,i\  {ud,i  -  Ug,i)  + 


Qd  =  Nd  \  mdK  +  27rrdKgNu  \^Td  -  Tgj 


Wd  =  Nd 


TT 


mdUd,iUd,i  +  ■^r^CDpg\ud,i  -  iud,i  -  Ug^i)  Ud,i  +  ^^Fd-^Ud,i 


(4.13a) 


(4.13b) 


(4.13c) 

(4.14a) 

(4.14b) 

(4.14c) 

(4.14d) 


where  Nd  is  the  number  density  of  the  Eulerian  dispersed  phase. 

The  source  terms  from  the  Lagrangian  and  Eulerian  dispersed  phases  provided  through 
Eqs.  (4.15)  and  (4.14)  can  be  combined  through  Eq.  (4.7)  to  yield  source  terms  for  the  gas 
phase.  The  governing  system  of  equations  given  by  Eqs.  (4.6),  (4.10)  and  (4.13),  for  the 
gas  phase,  Lagrangian  disperse  phase  and  Eulerian  disperse  phase,  respectively  are  complete 
once  the  subgrid-scale  models  are  specified,  which  are  briefly  described  in  the  next  section. 


4.2  Coupling  Between  the  Phases 

As  shown  in  Eq.  (4.7),  the  source  term  needed  for  the  gas  phase  transport  equations  given  by 
Eq.  (4.6)  comprises  of  source  terms  corresponding  to  the  Lagrangian  and  Eulerian  dispersed 
phases.  The  source  term  contribution  from  the  Lagrangian  dispersed  phase  are  obtained 
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through 


Pv  ~  AT/  ‘^p,n'fPp,n, 


1 


n=l 

N 


n=l 


rrip^nUp^^^n  + 


7T 


N 

'\ 

n. 


^g,i,n\  {}^p,i,n  ^g,i,n 


,2,77-) 


n=l 


^p^nhy^n  [Tp^n  ^g,n 


N 


^p,r 


77-=l 


^4  3  dpg^^ 

'l^p,nUp,i,nUp,i,n  +  ^^'''p,n  q^,  '^p,i,n 


^rri.n. 


/I  2  /^  I  ~  ff  \  f 

4”^^p,n^-D,nP3,n |^p,i,n  Pg,i,n  ^g,i,n\  \^P,i,i^  '^g,i,nj  '^p,t,n 

where  hy  is  the  enthalpy  change  associated  with  the  mass  transfer. 


(4.15a) 

(4.15b) 

(4.15c) 

(4.15d) 

(4.15e) 

(4.15f) 


4.3  Subgrid  Terms  and  Closure 

The  density-weighted  filtering  of  the  governing  equations  leads  to  appearance  of  unclosed 
terms,  also  referred  to  as  the  subgrid-scale  terms,  which  require  further  closure  approxima¬ 
tions  to  obtain  a  closed  system  of  equations.  As  mentioned  before,  these  terms  include  the 
subgrid  stress  tensor  subgrid  enthalpy  flux  subgrid  viscous  work  a®^®,  subgrid 

convective  species  flux  subgrid  diffusive  species  flux  d®®/^  and  the  subgrid  heat  flux 

%^ik  [54,  55].  Some  of  these  terms  can  be  closed  using  the  models  developed  for  non-reacting 
or  reacting  gas  phase  problems.  Here,  we  briefly  describe  some  well  known  models  that  are 
used  for  closure  of  subgrid-scale  terms. 

The  subgrid  stress  and  heat  flux  terms  are  typically  closed  following  the  Boussinesq 
approximation,  where  an  eddy  viscosity  type  gradient  closure  is  employed.  Such  type  of  clo¬ 
sure  is  very  popular  and  a  model  is  required  for  the  eddy  viscosity.  Both  algebraic/dynamic 
Smagorinsky  (ASM/DSM)  model  [56,  57]  and  the  model  for  subgrid  kinetic  energy  [58,  59] 
are  very  popular.  For  example,  when  using  the  subgrid  kinetic  energy  closure,  a  trans¬ 
port  equation  for  needs  to  be  solved,  which  is  given  by 

d _  d  , _ _  .  _  d  ( _ dk^^^  \  _ 

—ag  Pgk  +  —  [ag  Pg  Ug^ik  ^  )  -  agPksgs  +  —  \  agPgUt-g^  j  -  agDksgs  +  F^  k- 

(4.16) 

Here,  Pksgs  and  Dk^gs  denote  respectively,  the  production  and  dissipation  of  which  can 
be  obtained  through 


p  _  ^sgs^^g.i 

dxj  ’ 


Dksgs  —  C^p^ 


g^sgs-^ 


1.5 


(4.17) 
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where  Ut  is  the  subgrid  eddy  viscosity,  and  is  modeled  as  Ut  =  CuVk^A.,  where  A  is  the 
local  hlter  width.  The  constants  Ciy  and  Q  are  obtained  theoretically  as  0.067  and  0.916, 
respectively  [54]  but  can  be  computed  using  a  dynamic  procedure,  referred  to  as  the  locally 
dynamic  kinetic  energy  model  (LDKM)  [60,  55]  that  can  be  used  locally  without  requiring 
any  ad  hoc  averaging.  With  the  eddy  viscosity  model,  the  subgrid  stress  tensor  is  closed 
through 


rr-SQS 

9,V 


1  \  2_ 

—  2PgZ/f  I  Sg^ij  —  —Sg^hk^ij  j  +  3^9^  ^ 


(4.18) 


where  Sg^ij  is  the  rate  of  strain  tensor  given  by 

(4. 

—rsgs 

In  Eq.  (4.16)  ^  is  the  source  term  due  to  the  Lagrangian  and  Eulerian  dispersed  phases 

that  can  be  closed  exactly  [61,  62]. 

The  subgrid  total  enthalpy,  can  also  modeled  using  the  eddy  viscosity  and  gradient 

diffusion  assumption  through 


~  2 


1  /  du 


9,1- 


du. 


'90 


dxj 


dxj 


Tjsgs  _ 
^90 


-  9Hg 

^3  pr^  Qx-  ’ 


(4.20) 


where  Hg  is  the  filtered  total  enthalpy,  given  by  Hg  =  hg  +  \ug^iUg^i  +  where  hg  = 
J2k=i  hg,kYg,k-  The  turbulent  Prandtl  number  Prt  can  be  dynamically  computed,  but  is 
typically  assumed  to  be  unity.  In  a  similar  manner,  the  subgrid  convective  species  flux, 
is  modeled  using  the  gradient  diffusion  assumption  through 


ysgs  _  Pg^t  ^Yg^k 

Set  dxi  ’ 

where  Sct  is  the  turbulent  Schmidt  number.  It  can  also  be  dynamically  computed;  however, 
it  is  typically  assumed  to  be  unity  in  many  studies. 

The  other  subgrid  terms,  i.e.,  O^Pf,  and  are  neglected  in  most  studies  [54,  63].  Note 
that  in  the  Eq.  (4.18)-(4.21),  the  subscript  ‘p’  can  be  replaced  with  ‘d’  to  obtain  appropriate 
closures  under  the  assumption  that  the  Eulerian  dispersed  phase  behaves  as  a  pseudo  gas 
phase.  The  so-called  pseudo-fluid  assumption  is  not  very  easy  to  justify  and  therefore,  the 
closures  for  dispersed  phase  especially,  in  the  context  of  dense  flows  is  unknown. 

A  closure  model  for  the  reaction  rate  is  also  required  for  the  combustion  problem.  Subgrid 
turbulence-chemistry  interaction  models  have  remained  a  challenging  task  even  for  gas  phase 
combustion  and  these  problems  are  made  even  more  challenging  when  spray  evaporation, 
mixing  and  combustion  have  to  be  included.  Due  to  the  discrete  nature  of  the  particles 
vaporizing  and  mixing  is  local  and  therefore,  combustion  can  occur  in  multitude  of  manner 
ranging  from  non-premixed  to  fully  premixed  type.  Partially  premixing  is  the  norm  in  spray 
combustion  system  and  therefore,  subgrid  closures  for  turbulence-chemistry  interaction  need 
to  address  the  complete  regime.  At  this  time  there  are  not  that  many  options  for  this  goal. 
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4.4  Hybrid  Solver  transition 

As  two  solvers  for  the  simulation  of  dispersed  phase  are  available  in  LESLIE,  a  hybrid 
approach  has  been  developed  as  follows:  the  dense  phase  is  simulated  using  the  Eulerian 
particle  phase  approach,  while  the  diluted  phase  is  simulated  using  the  Lagrangian  phase. 

The  computations  using  Eulerian  particle  phase  are  fast  and  are  accurate  for  dense  to 
marginally  dense  flows.  Lagrangian  particle  results  are  accurate  over  a  wide  range  of  dis¬ 
persed  phase  volume  fractions  and  particle  size  distributions  and  this  technique  can  easily 
handle  polydisperse  flows.  Moreover,  as  every  particle  can  be  tracked  accurately  and  effi¬ 
ciently,  the  combustion  and  evaporation  of  Lagrangian  particles  can  be  studied  more  easily. 
The  Lagrangian  solver  with  dense  correction  can  also  be  used  for  heavily  loaded  particle 
flows,  but  it  can  be  very  expensive  when  the  number  of  particles  increases.  In  order  to  take 
advantage  of  both  methods,  a  hybrid  method  has  been  developed.  It  uses  the  Eulerian  ap¬ 
proach  in  dense  regions  where  the  liquid  volume  fraction  is  high,  and  Lagrangian  approach 
when  the  loading  is  much  lower.  The  technique  to  transfer  information  from  the  Eulerian 
dispersed  phase  to  Lagrangian  phase  is  determined  by  three  parameters,  namely,  the  transfer 
volume  fraction  ctt,  the  transfer  fraction  /t  ,  and  the  number  of  particles  per  parcel,  Pp 
.  The  transfer  volume  fraction  is  the  threshold  value  to  switch  from  the  Eulerian  solver  to 
the  Lagrangian  one  {a  <  ar)-  It  indicates  the  region  where  the  transition  from  EE  to  EL  is 
desired.  If  this  parameter  is  very  low,  the  solver  remains  purely  Eulerian  and  no  Lagrangian 
particles  are  created.  On  the  other  hand,  if  alphax  is  set  to  a  relatively  high  value,  all  the 
Eulerian  particles  are  transferred  into  the  Lagrangian  solver.  It  needs  to  be  between  0  and 
1.  The  choice  of  this  parameter  can  be  tricky,  however,  it  depends  on  the  computational 
cost  and  the  required  accuracy. 

The  transfer  fraction  is  the  fraction  of  particles  in  a  computational  cell  transferred  from 
Eulerian  dispersed  phase  to  Lagrangian  dispersed  phase.  When  the  particles  are  added  to 
the  Lagrangian  dispersed  phase,  the  volume  fraction  of  the  Lagrangian  dispersed  phase,  ftp, 
is  incremented  based  on  the  number  of  particles  added.  The  number  of  particles  added  to  the 
Lagrangian  dispersed  phase  is  given  as:  Np  =  As  the  particles  are  removed  from 

the  Eulerian  dispersed  phase,  the  Eulerian  liquid  volume  fraction  is  recomputed  to  conserve 
mass. 

The  mass,  the  momentum,  the  energy,  the  velocity,  the  temperature,  the  radius  and  the 
position  of  each  particle  transferred  to  the  Lagrangian  dispersed  phase  are  set  based  on  the 
values  from  the  Eulerian  dispersed  phase  in  a  given  computational  cell.  If  /t  =  1,  then  all 
the  particles  in  a  computational  cell  are  transferred  and  if  /^  =  0  no  particle  is  transferred. 
A  value  between  0  and  1  is  chosen  to  ensure  stability  instead  of  transferring  all  the  particles. 
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CHAPTER  V 


RESULTS  AND  DISCUSSIONS 


5.1  Validation 

5.1.1  Sedov  test  cases 

To  assess  the  capability  of  the  code  to  accurately  solve  blast  waves,  a  first  set  of  simulations 
has  been  performed  and  compared  to  the  Sedov  analytical  solutions.  These  simulations 
assume  a  point  explosion  with  a  pressure  field  of  10“^  Pa.  A  high  pressure  blob  is  initialized 
at  the  center  of  the  domain  with  a  1000  Pa  value  in  a  5  mm  radius  disk.  The  density  is 
equal  to  1  kg/m^  everywhere.  Outflow  conditions  are  imposed  at  the  boundaries  to  avoid 
reflections.  The  simulated  physical  time  is  100  ms.  The  evolution  of  the  pressure  peak  and 
shock  radius  is  compared  on  three  different  uniform  grids  (256  x  256,  512  x  512,  1024  x 
1024),  and  in  3  dimensions,  inside  a  3D  sector  on  three  different  meshes  (250  x  45  x  45,  500 
X  45  X  45,  1000  X  45  X  45),  where  the  two  angles  (p  and  6  are  45  degrees. 

The  evolution  of  the  shock  radius  and  the  pressure  peak  is  presented  in  Fig.  5.1.  The 
comparison  is  done  using  all  the  uniform  meshes  presented  in  this  section  with  an  additional 
simulation  computed  using  AMR  (Adaptative  Refinement  Method)  around  the  shock.  An 
empirical  law  is  compared  to  these  profiles.  Results  give  good  agreements  for  both  the  shock 
radius  and  the  pressure  peak  evolution. 

The  quantitative  comparison  is  provided  in  Fig.  5.2,  comparing  the  Shock  radius  and 
the  pressure  peak  evolution  with  time.  These  profiles  are  compared  to  an  extra  simulation 
using  AMR  around  the  shock  radius,  and  an  empirical  law  from  Sedov  derivations.  From 
this  study,  a  very  good  agreement  is  found. 

5.1.2  Boiko  test  case 

In  this  test,  the  EL  solver  is  validated  by  simulating  the  interaction  of  a  marginally  dense 
particle  cloud  with  a  shock  wave.  The  setup  for  this  test  is  based  on  the  experiments 
performed  by  Boiko  et  al.  [1].  A  shock  tube  of  length  6.5  m  and  cross  section  52  mm  x 
52  mm  was  used  in  the  experiment  and  an  acrylic  plastic  particle  cloud  of  initial  volume 
fraction  3.0%  was  impacted  by  a  Mach  2.8  shock.  The  radius  of  each  particle  is  150  jim 
and  the  density  is  1200  kg/m^.  For  simulation,  a  domain  with  dimensions  same  as  that  of 
the  shock  tube  is  considered  and  is  resolved  using  as  grid  of  size  1250  x  10.  The  particle 
cloud  is  placed  at  a  location  2  m  from  the  high  pressure  end  of  shock  tube.  The  initial 
width  of  the  particle  cloud  is  13  mm  and  the  cloud  has  74600  particles.  The  driver  section 
is  filled  with  He  and  is  initialized  with  high  pressure  p4  obtained  from  Eqn.  (5.1)  such  that 
a  Mach  2.8  shock  interacts  with  the  particle  cloud.  The  driven  section  is  filled  with  air  at 
10^  Pa.  The  driven  section,  the  driver  section  and  the  particle  cloud  are  initially  at  298K. 
All  boundaries  of  the  shock  tube  are  set  to  be  no-flux  boundaries  except  the  outflow  where 
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(a)  Shock  Radius 


(b)  Pressure  peak  evolution 

Figure  5.1:  Shock  radius  and  pressure  peak  evolution  with  time  for  2D  simulations:  all  cases 
are  compared  with  empirical  law  and  an  extra  simulation  done  using  Adaptative  Refinement 
method 
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Figure  5.2:  Shock  radius  and  pressure  peak  evolution  with  time  for  3D  simulations:  all  cases 
are  compared  with  empirical  law  and  an  extra  simulation  done  using  Adaptative  Refinement 
method 
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Figure  5.3:  Dispersion  of  particle  cloud  by  a  shock  wave.  The  experimental  data  is  from 
[1] 


the  flow  properties  are  extrapolated.  The  acceleration  of  the  particle  is  computed  using  the 
drag  law  provided  by  Boiko  et  al.  [1], 


JU  ^  27iMf  -  (71  -  1)  r  _  (74  -  1)  ai  _  J_\ 

Pi  7i  + 1  L  (71  + 1)  04  V  ^  MiJ 

The  dispersion  of  the  particle  cloud  after  interaction  with  the  shock  wave  is  shown  in 
Fig.  5.3.  Here,  unlike  the  single  particle  dispersion  case,  DEM  is  turned  on  and  hence  the 
gas-phase  fluxes  are  influenced  by  the  volume  of  the  particles  in  the  flow.  The  case  shows 
the  accuracy  of  the  solver  in  simulating  the  particle  dispersion  for  marginally  dense  clouds. 

The  Boiko  test  is  repeated  with  EE  method  using  a  grid  with  resolution,  A  =  250.0  /rm 
and  the  variation  of  the  location  of  the  particle  cloud  with  time  is  compared  with  the  results 
from  EL  method.  While  EL  method  accurately  captures  the  particle  cloud  motion,  EE 
method  shows  an  agreement  in  initial  stages  of  the  cloud  dispersal.  As  the  volume  fraction 
of  the  particle  cloud  decrease  from  0.03  to  below  the  dilute  limit  (i.e.  0.01),  the  error  in  the 
particle  cloud  location  obtained  using  EE  scheme  increases  from  1.5  %  to  33.0  %.  Since, 
the  particle  cloud  location  is  dependent  on  the  accurate  estimation  of  the  cloud  extremities, 
the  error  with  EE  increases  in  this  case.  Note  that  the  dilute  dispersion  with  EE  is  not  as 
accurate  as  EL  as  EL  accounts  for  the  motion  of  each  particle  and  accurately  determines  the 
particle  cloud  edges. 
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Figure  5.4:  Comparision  of  pressures  upstream  (US)  and  downstream  (DS)  of  a  dense 
particle  cloud.  The  experimental  data  is  from  [2],  The  results  with  EE,  EL  and  EE-EL  are 
shown.  Also,  EL  case  without  DEM  is  shown. 


5.1.3  Rouge  test  case 

The  propagation  of  a  shock  wave  through  a  dense  particle  cloud  is  simulated  using  both  EE 
and  EL  methods  to  validate  both  the  methods  and  investigate  the  importance  of  the  gas- 
phase  flux  correction  based  on  the  particle  volume  fraction.  The  setup  for  this  test  is  based  on 
the  experimental  conflguration  described  by  Rouge  et  al.  [2] .  A  vertical  shock  tube  of  length 
6  m  and  cross  section  13  cm  x  13  cm  has  been  employed  in  the  experimental  investigations. 
Glass  particles  of  radius  750.0  fim  are  placed  in  the  shock  tube  and  are  initially  supported 
by  a  membrane.  A  Mach  1.3  shock  is  allowed  to  interact  with  the  particles  which  initially 
form  a  bed  of  thickness  2  cm  and  volume  fraction  65%.  For  these  simulations,  the  geometry 
of  the  domain  is  set  to  be  same  as  the  experimental  conflguration  and  is  resolved  using  a  grid 
of  size  600  x  13  x  13.  In  order  to  perform  simulations  using  EL  solver,  15540  computational 
particles,  i.e.,  parcels,  are  initialized  with  8  particles  per  parcel.  The  shock  tube  is  filled 
with  air  and  the  pressure  in  the  driver  section  is  set  based  on  Eqn.  (5.1).  The  pressure 
in  the  driven  section  is  10^  Pa,  initially.  The  initial  temperature  of  the  particles  and  the 
gas  is  at  298  K.  The  velocity  of  the  particles  is  computed  based  on  the  quasi-steady  drag 
relation  described  by  Crowe  et  al.  [64].  All  walls  of  the  shock  tube  are  set  to  be  no- flux 
boundaries  except  the  outflow  (the  end  away  from  the  high  pressure  driver  section)  where 
the  flow  properties  are  extrapolated. 
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The  pressure  at  locations  11  cm  below  the  particle  bed  and  72  cm  above  the  particle  bed 
are  recorded  and  compared  with  the  experimental  results.  As  the  shock  interacts  with  the 
dense  particle  bed,  a  transmitted  shock  propagates  through  the  particle  cloud  and  a  reflected 
shock  wave  propagates  upstream.  In  order  to  capture  the  pressure  accurately,  the  gas-phase 
flux  correction  based  on  the  volume  fraction  of  the  dispersed  phase  is  very  important.  This 
is  demonstrated  in  Fig.  5.4.  In  the  case  without  DEM,  the  upstream  pressure  is  under 
predicted  and  the  downstream  pressure  is  over  predicted.  Also,  the  downstream  shock  wave 
propagates  through  the  particle  cloud  at  relatively  greater  speed  and  reaches  the  upstream 
pressure  trace  point  nearly  0.4  ms  earlier.  Thus,  DEM  is  important  to  predict  the  gas- 
phase  and  the  dispersed  phase  properties  and  its  significance  is  increased  with  increase  in 
the  dispersed  phase  volume  fraction. 

Eor  this  test,  the  numerical  simulations  are  performed  using  both  EE  and  EL  methods. 
The  cases  with  EE  are  performed  with  and  without  the  contribution  from  granular  friction 
{pf)  and  granular  viscous  dissipation  terms  (7  and  (f)).  Eigure  5.4  shows  that  the  effect  of 
friction  and  the  granular  viscous  terms  is  not  significant.  Both  EE  and  EL  methods  produce 
similar  results  and  are  in  good  agreement  with  the  experimental  results.  The  results  with 
the  combined  EE-EL  method,  where  the  dispersed  phase  is  transitioned  from  EE  method  to 
EL  method  if  the  Eulerian  dispersed  phase  volume  fraction  reaches  0.01,  are  also  in  good 
agreement  with  the  results  from  the  experiment  and  the  simulations  using  pure  EE  and  pure 
EL  methods. 


5.2  Blast  Wave  in  Room  with  Vent 

5.2.1  Simulation  setup 

In  this  study,  we  simulate  the  propagation  of  a  blast  wave  from  a  spherical  charge  in  an 
enclosure  with  a  vent.  The  explosive  is  initially  set  in  the  middle  of  the  enclosure.  A  two- 
dimensional  representation  of  the  setup  is  simulated.  The  case  is  initialized  with  detonation 
wave  profiles  for  the  pressure,  velocity,  density,  and  temperature  which  are  obtained  using 
the  Gas-Interpolated-Solid  Stewart-Prasad-Asay  (GISPA)  method  for  the  condensed  phase 
detonation  process.  The  energetic  material  in  consideration  is  TNT  and  the  charge  is  11.63 
cm  in  diameters.  The  room  measures  20  charge  diameters  in  length,  16  diameters  in  height, 
and  the  vent  has  a  width  of  2  diameters  located  in  the  center  of  the  ceiling  (top  boundary) ,  as 
shown  in  the  schematic  diagram  in  figure  5.5(a).  The  simulation  captures  the  pressure  front 
and  subsequent  reflections  from  the  walls.  These  develop  locations  with  wave  focusing  where 
pressure  spikes  can  be  observed.  The  simulation  makes  use  of  adaptive  mesh  refinement 
(AMR)  to  better  resolve  the  thin  blast  front. 

5.2.2  Results  and  discussions 

Eigure  5.6  plots  the  time  evolution  of  the  pressure  at  the  location  of  the  probes  as  described 
in  5.5(b).  Eigures  5.6  a,  b,  and  c  compare  corresponding  probes  mirrored  with  a  line  crossing 
horizontally  in  the  center  of  the  domain.  Up  until  0.4ms,  the  probes  show  excellent  agreement 
demonstrating  a  symmetry  in  the  flow  behavior.  This  is  expected,  since  at  this  stage,  the 
expanding  blast  wave  is  symmetric,  and  has  yet  to  reach  the  vent,  which  is  the  source  of 
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20  CD 


(a)  Schematic  diagram 


Figure  5.5:  Schematic  diagram  of  domain  and  location  of  sensors. 
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asymmetry  in  the  geometry  The  three  sensor  show  a  similar  behavior  in  the  number  of 
peaks  depicting  the  wave  front  passage  and  subsequent  reflections.  Figure  5.6  d  compares 
the  pressure  traces  along  the  vertical  centerline.  The  plot  depicts  an  interesting  pressure 
focusing  in  the  center  of  the  domain  with  higher  peaks  than  its  surrounding. 

The  propagation  of  the  blast  wave  undergoes  multiple  characteristics  that  can  be  split 
here  into  three  stages.  The  first  stage  can  be  seen  as  the  outward  propagation  of  the  blast 
wave  from  the  center  charge.  The  blast  wave  maintains  the  initial  spherical  shape,  which  is 
imprinted  from  the  initial  charge  shape.  This  can  be  seen  in  parts  a  and  b  of  figures  5.7  and 
5.8.  The  second  stage  appears  as  the  waves  get  reflected  from  the  walls.  The  wave  reaches 
the  top  and  bottom  walls  at  the  center  and  a  pressure  peak  propagates  along  the  walls.  The 
reflected  waves  then  coalesce  at  the  corners  propagating  inwards.  The  propagation  now  takes 
the  imprinted  shape  of  the  enclosure.  The  rectangular  shape  can  be  most  clearly  seen  in 
part  d  of  figures  5.7  and  5.8.  However,  the  rectangular  shape  is  not  complete,  particularly  at 
the  top  part,  where  the  vent  imprint  is  noticeable.  The  vent  allows  pressure  relief  and  thus, 
a  delay  in  the  wave  propagation  in  this  region  is  observed.  The  third  stage,  seen  in  parts  e 
and  f  of  figures  5.7  and  5.8  show  a  merging  of  the  reflected  waves  that  are  now  propagating 
back  towards  the  walls.  This  wave  propagation  continues  while  the  waves  strength  reduces 
as  shows  in  the  time  traces  in  figure  5.6. 

5.3  Blast  wave  in  a  Room  Configuration 

5.3.1  Simulation  Setup 

In  this  case,  the  computational  domain  is  a  rectangular  shape  single  room  which  measures 
3[m]  in  width  and  depth,  and  l[m]  in  height.  The  imposed  boundary  condition  is  slip 
walls  everywhere.  The  charge  is  cylindrical  in  shape,  measuring  0.02286[m]  radius  and 
0.2286[m]  length,  and  placed  at  the  center  of  computational  domain.  Two  different  methods 
of  blast  wave  initialization  are  applied  to  the  simulation  in  order  to  validate  them.  The  first 
method  is  a  full  burnt  simulation,  which  imports  one-dimensional  detonation  profiles.  Thus, 
the  physical  quantities,  such  as  pressure,  temperature,  etc,  are  not  distributed  in  the  axial 
direction  of  charge.  On  the  other  hand,  the  second  method,  uses  a  center  burn  initialization 
where  the  detonation  propagates  from  the  center  to  the  extremities  of  charge.  This  case  is 
expected  more  faithfully  represent  the  real  life  burning  compared  to  the  full  burnt  case.  For 
full  burnt  case,  the  charge  is  composed  by  TNT,  while  HMX  is  used  for  center  burnt  case. 

To  measure  the  influence  of  Particles,  the  same  computational  grid  is  utilized.  The  size 
of  explosive  charge  is  kept  the  same,  and  HMX  with  center  burnt  initialization  is  applied. 
The  particles  are  initially  placed  at  charge,  and  propagates  into  the  room  driven  by  the  blast 
wave.  The  parameters  used  for  particle  initialization  is  shown  in  Table  5.1. 

5.3.2  Results  and  Discussion 

At  first,  the  grid  resolution  is  investigated  to  guarantee  the  accuracy  of  resolution.  The 
coarse,  middle,  and  fine  grid  which  respectively  has  about  20,  40,  and  80  nodes  are  used  for 
this  verification.  The  center  burnt  case  is  used  for  initialization.  The  pressure  history  at 
three  different  locations  are  shown  in  Figure  5.9.  The  coarse  grid  shows  a  small  deviation 


23 


(a)  (b) 


time  (ms)  time  (ms) 

(C)  (d) 

Figure  5.6:  Time  evolution  of  pressure  at  the  different  locations. 


Table  5.1:  Initial  setup  for  particles. 


Case 

Number  of 
Particles 

Radius  [jim] 

Number  of 
particles  per 
parcel 

Volume  fraction 

A 

0 

- 

- 

- 

B 

100,000 

2.0 

1 

8.93e-2 

C 

100,000 

20.0 

1 

8.893e-6 

D 

100,000,000 

2.0 

1 

8.93e-6 

E 

100,000,000 

20.0 

1 

8.93e-3 
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Figure  5.7:  Pressure  contour  plots. 
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Pressure  (1)  Pressure  (1)  Pressure  (1) 


Figure  5.8:  Velocity  contour  plots. 
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Figure  5.9:  Pressure  history  at  three  different  locations.  All  locations  are  in  center  plane 
which  is  parallel  to  the  floor.  The  distance  from  center  axis  of  charge  d  varies;  (a)  d  = 
0.02286[m],  (b)  d  =  0.5[m],  and  (c)  d  =  1.0[m]. 


from  the  other  two  cases  althongh  the  same  tendency  is  captnred.  The  middle  and  fine  case 
shows  good  agreement  in  qnantitative  sense,  which  indicates  that  the  grid  is  fine  enough  in 
middle  size  grid.  Therefore,  the  middle  size  grid  is  used  for  all  following  studies. 

Figure  5.10  to  5.13  show  the  evolution  of  the  blast  wave  for  CASE  A.  Initially  the  center 
part  of  charge  is  burnt  and  has  high  pressure,  whereas  the  other  part  of  charge  is  unburnt 
solid  phase.  Therefore,  the  jet  going  to  outside  is  formed  at  first.  The  solid  phase  prevent  the 
shock  wave  from  propagating  into  the  vertical  direction,  which  resnlts  in  high  pressure  region 
at  the  contact  surface  of  charge  and  detonation  products.  As  time  passes,  the  detonation 
propagates  in  the  vertical  directions.  Hence,  the  pressure  distributes  along  the  vertical 
directions,  and  has  a  curvature.  The  combustion  of  the  charge  completes  around  t  =  14[s]. 
After  the  combustion,  the  detonation  products  start  expand  to  the  vertical  direction  with 
significant  velocity.  This  shock  wave  expanding  to  the  vertical  direction  initially  has  a  flat 
shape  due  to  the  influence  of  charge  shape.  Also,  there  is  a  slight  discontinnity  between 
shock  expanding  to  the  horizontal  and  vertical  directions,  which  is  similar  to  the  one  in  full 
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burnt  case. 

For  the  particle  case,  the  evolutions  of  blast  wave  shape  for  case  A  to  E  are  shown  in 
Fignre  5.14.  In  these  figures,  there  are  some  perturbations  at  the  top  and  side  of  the  shock. 
Althongh  these  fignres  show  that  lighter  particles  tends  to  move  faster  than  heavy  one,  there 
is  no  significant  difference  between  all  cases  in  quantitative  sense.  Also,  the  pressure  histories 
seem  identical  for  all  cases,  which  are  shown  in  Figure  5.15.  This  indicates  that  the  numerical 
diffusion  due  to  the  numerical  scheme  or  grid  resolution  is  so  large  that  the  perturbations 
indnced  by  particles  are  dissipated. 


5.4  Blast  Wave  in  Rooms  with  Hallways 

5.4.1  Simulation  Setup 

The  initial  detonation  profiles  were  obtained  from  Dr.  Douglas  Nance  (Eglin  Air  Force  Base) 

[65] .  This  is  based  on  a  one-dimensional  Direct  Nnmerical  Simulation  (DNS)  employing  the 
Gas-Interpolated-Solid  Stewart-Prasad-Asay  (GISPA)  method  for  the  detonation  process 

[66]  .  This  method  permits  time-accurate  simulation  of  detonation  from  the  time  of  the  initial 
shock  through  the  completion  of  the  explosive  burn.  The  robustness  of  the  GISPA  algorithm 
is  emphasized  by  its  ability  to  capture  the  reaction  zone  as  well  as  the  Von  Neumann  spike. 
The  GISPA  method  is  based  upon  the  reactive  Euler  equations  [67]. 

f  +  f  =  +  (5.2) 

Here  U  =  (p,  u,  E,  pX)^  is  the  vector  of  conserved  variables.  A  denotes  the  reaction 
progress  state  variable,  and  F  =  {pu,  pu^+p,  u(pE+p),  puX)'^)  is  the  fiux  vector.  The  source 
term  So  mathematically  corrects  the  one-dimensional  equations  for  non-planar  coordinate 
systems. 


Sg  ^  -^{pu,  pu^  +  p,u{pE  +  p),  puX)'^)  (5.3) 

Here  x  denotes  distance,  and  j  is  set  to  0  for  planar,  1  for  cylindrical,  and  2  for  spherical. 
The  progress  of  the  detonation  is  governed  by  a  reaction  rate  expression  Sr^,  which  can 
take  different  forms  for  different  explosives  [68].  These  equations  are  solved  with  the  use  of 
appropriate  equations  of  state  for  both  the  condensed  explosive  and  the  detonation  products. 
For  the  condensed  explosive,  the  Hayes  equation  of  state  is  used  [69] ,  while  the  JWL  equation 
of  state  is  used  for  the  detonation  products.  The  Glaister’s  [70]  version  of  the  Roe  scheme  is 
used  with  MUSCL  reconstruction  for  solving  the  equations  to  obtain  the  initial  detonation 
profile.  The  detonation  initialization  based  on  the  GISPA  method  is  more  realistic  than  other 
ways  of  initialization,  snch  as  the  “programmed  burn”  algorithm  [71],  and  the  constant 
volume  explosion  initialization.  In  the  programmed  burn  algorithm,  the  detonation  wave 
speed  needs  to  be  known  a  priori,  which  may  not  always  be  the  case.  In  the  constant  volume 
explosion  initialization,  since  the  pressnre  field  is  assnmed  constant,  the  early  momentum 
transfer  characteristics  from  the  gas  to  the  particles  can  be  erroneous.  Since  the  GISPA 
algorithm  is  based  on  first  principles,  we  believe  that  this  initialization  is  more  realistic  than 
the  other  procedures. 
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(e)  t  =  16/is  (f)  t  =  16/is 


(g)  t  =  20fis 


(h)  t  =  20/is 


Figure  5.10:  Blast  wave  formation  of  full  burnt  case.  Left:  pressure,  right:  velocity. 
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(a)  t  =  0/i5 


(b)  t  =  Ojis 


(c)  t  =  8/is  (d)  t  =  8/is 


(e)  t  =  16/i5 


(f)  t  =  16/i5 


(g)  t  =  20/us 


(h)  i  =  20/us 


Figure  5.11:  Blast  wave  formation  of  full  burnt  case.  Left:  Density,  right:  Mass  Fraction 
of  CO. 
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(a)  t  =  0/is  (b)  t  =  0/is 


(e)  t  =  16/is 


(f)  t  =  16/is 


(g)  t  =  20/is  (h)  t  =  20/is 


Figure  5.12:  Blast  wave  formation  of  center  burnt  case.  Left:  pressure,  right:  velocity. 
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(a)  t  =  Qjis  (b)  t  =  0/is 


(c)  t  =  8/i5 


(d)  t  =  Sjis 


(e)  t  =  IQfis 


(f)  t  =  16/i5 


(g)  t  =  20/us 


(h)  t  =  20/us 


Figure  5.13:  Blast  wave  formation  of  center  burnt  case.  Left:  density,  right:  mass  fraction 
of  DHMX. 
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♦  # 


(a)  Case  A 

(b)  Case  B 

♦  ♦  ♦ 

(c)  Case  C 

*  ♦  # 

(d)  Case  D 

*  # 

(e)  Case  E 

Figure  5.14:  Evolution  of  blast  wave.  Contour  corresponds  iso-surface  of  Ydhmx  =  0.5 
and  particles  distributions  are  also  shown.  From  left  to  right,  hgure  corresponds  t  =  8,  16 
20,  and  30  [s] 
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(c) 


Figure  5.15:  Pressure  profiles  for  case  A  to  E.  All  locations  are  in  center  plane  which  is 
parallel  to  the  floor.  The  distance  from  center  axis  of  charge  d  varies;  (a)  d  =  0.02286[m], 
(b)  d  =  0.5[m],  and  (c)  d  =  1.0[m]. 
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Table  5.2:  Soot  burning  mechanism  implemented  in  current  study 
Number  Elementary  reactions 

1  C{S)  +  O.5O2  =  CO 

2  C0  +  O.5O2  =  CO2 


Table  5.3:  Simplified  reaction  mechanism  for  C4 


Number 

Elementary  reactions 

1 

OH  +  0  =  H  +  O2 

2 

OH  +  H  =  H2  +  0 

3 

OH  +  OH  =  H20  +  0 

4 

OH +  H2  =  H20  +  H 

5 

H  +  H  +  M  =  H2  +  M 

6 

H  +  0  +  M  =  OH  +  M 

7 

0  +  0  +  M  =  02  +  M 

8 

H  +  OH +  M  =  H20  +  M 

9 

OH  +  CO  =  C02  +  H 

10 

C0  +  0  +  M  =  C02  +  M 

5.4.2  Chemical  Kinetics 

MRTF  tests  are  conducted  with  C4  explosive.  Additionally  in  this  work  NM  is  used  to  carry 
out  comparative  study  to  check  the  charge  effect  on  pressure  signals.  C4  is  composed  of  91% 
RDX  (C^HeNeOe)  and  9%  binders  such  as  DOR  {H3sC2^04),  PIB  (//gQ),  and  fuel  oil.  The 
constituents  of  the  fireball  for  C4  and  NM  charge  are  defined  as: 

C4  Detonation  :  2.713H20  +  1.0732H2  +  2.608N2  +  2.608CO  +  1.252C(S) 

NM  Detonation  :  0.2296N2  +  0.295H2O  +  0.459CO  +  0.0164H2 

Chemical  kinetics  mechanism  for  hydrocarbons  is  typically  very  complex  involving  nu¬ 
merous  species  reacting  with  each  other  intermittently  in  several  different  steps.  Currently, 
a  simple  two  step  mechanism  has  been  employed  in  the  study.  This  chemistry  assumes  soot 
to  be  present  in  the  gaseous  phase. 

However,  this  approach  has  its  limitation  and  therefore  an  implementation  of  a  more 
complex  and  detailed  chemistry  is  under  progress.  The  proposed  chemistry  mechanism 
involves  9  species  and  10  reversible  reactions  [72]. 


Table  5.4:  Simplified  soot  oxidation  path 
Number  Elementary  reactions 

1  C{S)  +  O.5O2  =  CO 

2  CO  +  O.5O2  =  CO2 
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Figure  5.16:  MRTF  test  structure  and  QSDD  instrumentation 


5.4.3  Experimental  and  Numerical  Configuration 

The  MRTF  or  the  multi  room  test  facility  is  a  five  room  concrete  structure,  as  shown 
in  Figure  5.40  [72],  There  is  an  internal  hall- way  which  connects  these  rooms.  Pressure 
and  temperature  gauges  are  installed  in  the  entire  structure  to  measure  various  explosive 
parameters. 

Room  3  is  the  cynosure  of  this  study  as  this  room  serves  as  the  static  detonation  room. 
The  room  is  fit  with  9  pressure  gauges  to  measure  the  shock  pressure  and  the  quasi-static 
pressure.  The  room  also  includes  3  temperature  gauges.  The  room  dimensions  are  shown 
in  Figure  5.17.  Its  a  21  ft.  X  19  ft.  X  9.5  ft.  room  with  a  10  vent  in  the  ceiling  at  the 
center.  There  is  a  4  ft.  X  4  ft.  door  which  connects  it  to  the  hallway.  Another  configuration 
including  the  hall-way  has  been  tested  to  make  sure  that  the  outflow  boundary  condition  at 
the  exit  of  the  room  is  not  making  any  difference  to  the  chamber  pressure  signals.  A  detail 
of  this  configuration  is  given  in  Figure  5.17.  Twenty  lbs  of  C4  is  detonated  at  the  center  of 
this  room.  The  uncased  cylindrical  charge  is  4.5  in  diameter  and  22.5  long. 

The  computational  grid  required  for  the  simulation  has  been  created  using  ICEM-CFD 
meshing  software.  The  multi-block,  structured  mesh  is  composed  of  2632  blocks  enumerating 
22  million  grid  points.  Figure  5.18  shows  the  block  distribution  and  the  mesh.  The  blocking 
has  been  done  in  such  a  fashion  that  it  enables  to  have  a  cylindrical  grid  at  the  center  of  the 
room  to  represent  a  cylindrical  charge  accurately. 
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A. 


Figure  5.17:  2a 


Dimension  details  of  Room  3,  2b  :  Room  3  with  hallway 


Figure  5.18:  Block  distribution  and  overall  mesh  details  of  Room  3  grid 
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Figure  5.19:  Section  view  of  the  0-grid  (Room  3) 


Figure  5.20:  Room  3  grid  with  hall  way 


0-grids  are  used  in  the  cylindrical  region  section  while  careful  grid  stretching  has  been 
implemented  to  gradually  relax  the  grid  size  towards  the  walls  of  the  room.  This  strategy 
saves  grid  points  and  thereby  reduces  computational  effort.  Figure  5.19  shows  the  section 
view  of  grid  with  core  section  detailing  the  0-grid. 

Grid  coarsening  technique  has  been  implemented  for  the  room  with  hall-way  configura¬ 
tion.  Figure  5.20  shows  the  details  of  this  configuration  and  grid. 

Wall  boundary  condition  is  used  on  all  walls  whereas  supersonic  outflow  condition  is 
implemented  at  the  vent-hole  exit  and  the  room-door  exit  plane. 

5.4.4  Results  and  Discussions 

In  this  section,  results  for  various  simulations  have  been  reported.  Table  5.5  summarizes  the 
run  matrix. 

Run  1  and  2  evaluate  the  effect  of  the  exit  boundary  condition.  Run  2  domain  is  limited 
to  the  exit  door  of  Room  3.  The  objective  of  comparing  these  two  runs  is  to  check  whether, 
wave  reflection  from  the  door  is  significant  to  perturb  the  pressure  signals  as  recorded  by 
the  pressure  transducers  fit  in  Room  3.  Most  of  the  earlier  studies  dealing  with  blast  or 
explosion  performed  at  CCL  have  dealt  with  spherical  charge.  In  this  study,  the  interpolation 
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Table  5.5:  Simulation  Matrix 


Run 

Configuration 

Explosive 

Reactive 

Charge  Shape 

1 

Room  3  with  Hall 

C4 

No 

Cylindrical 

2 

Room  3 

C4 

No 

Cylindrical 

3 

Room  3 

C4 

No 

Spherical 

4 

Room  3 

NM 

No 

Cylindrical 

5 

Room  3 

C4 

Yes 

Cylindrical 

methodology  to  get  the  initial  detonation  profiles  has  been  extended  to  cylindrical  charges  to 
calculate  the  initialization.  A  comparison  study  between  a  spherical  and  cylindrical  charge 
has  been  carried  out  to  understand  the  effect  of  charge  shape  on  pressure  signals  and  flow 
field.  A  comparative  study  between  two  different  charges  e.g.  NM  and  C4  has  also  been 
reported.  While  the  first  four  runs  are  carried  out  in  non-reactive  setup,  the  fifth  run  is 
performed  with  chemistry  switched  on.  Te  simulation  runs  were  started  with  a  grid  with  16 
million  grid  points.  There  after  a  systematic  grid  refinement  is  done  throughout  the  domain 
to  achieve  a  grid  independent  solution.  The  final  grid  has  22  million  grid  points.  It  should 
be  noted  that  the  study  has  been  carried  out  for  Room  3  in  a  non-reactive  domain  with  the 
purpose  of  saving  computational  time.  Once,  the  grid  independency  has  been  achieved,  the 
same  mesh  resolution  has  been  applied  to  the  hall  way. 

One  of  the  primary  objectives  of  the  project  is  to  evaluate  soot  particle  dispersion  char¬ 
acteristics  from  the  Euler-Lagrangian  solution.  Soot  particles  are  typically  of  3-3.5  micron 
size.  To  have  considerable  soot  mass  loading,  it  would  require  tracking  billions  of  particles. 
Eulerian-Lagrangian  calculations  are  computationally  expensive,  so  any  strategy  to  reduce 
domain  dimension  and/or  grid  points  is  cherished.  The  pressure  pulse  results  for  Room  3 
configuration  and  Room  3  with  hallway  configuration  has  been  presented  in  Eigure  5.21. 
Again,  it  should  be  noted  that  this  results  are  for  non-reactive  runs  and  the  grid  density  for 
this  run  is  intermediate  (16  million).  From  the  results  it  can  be  seen  that  there  is  insignifi¬ 
cant  (i5%)  difference  in  pressure  pulse  till  55  millisecond.  After  wards,  there  is  considerable 
amount  of  difference  in  the  pressure  pulse  behavior.  Room  3  shows  higher  value  afterwards 
primarily  because  of  pressure  wave  reflection  from  the  door  exit. 

Figure  5.22  shows  the  instantaneous  flow  field  of  pressure.  The  pressure  waves  leave  the 
exit  door  of  Room  3  and  get  reflected  from  the  walls  of  the  hall.  Adding  the  other  rooms 
to  this  domain  would  therefore  be  expected  to  completely  mitigate  the  wave  reflection  and 
thereby  affecting  the  pressure  signals. 

Figure  5.23  shows  the  soot  mass  fraction  plots  at  two  different  time  instances. 

5.4.5  Non  reactive  C4  Charge 

Figure  5.24  shows  the  comparison  of  pressure  time  history  between  experimental  data  and 
non-reactive  C4  explosion  in  Room  3  configuration  for  non-reactive  case  (without  after  burn¬ 
ing).  From  the  chart  it  is  visible  that  there  exists  considerable  amount  of  delay  for  the  first 
and  subsequent  shock  arrival  time.  Also,  the  peak  pressure  is  less  for  the  simulation  by 
about  50  psi.  As  time  progresses,  the  delay  in  shock  arrival  time  also  increases.  However, 
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Figure  5.21:  Comparison  of  the  pressure  signals  between  Room  3  and  Room  3  with  hallway 
configurations 


Figure  5.22:  Room  3  with  hall  configuration:  Instantaneous  pressure  fields  at  15.5  and 
34.2  milliseconds 


the  number  of  peaks  obtained  from  the  simulation  is  quite  satisfactory  as  it  matches  the 
experimental  data  closely. 

Figure  5.25  shows  instantaneous  flow  fields  of  pressure  for  Room  3  configuration  on  a 
section  plane  that  is  at  the  center  of  the  room.  The  explosion  initiates  a  pressure  wave 
which  hits  the  wall,  gets  reflected  and  moves  towards  the  center  of  the  room  colliding  with 
the  next  wave  generated  by  the  blast.  Every  pressure  peak  in  Figure  5.25  corresponds 
to  the  situation  when  two  different  pressure  waves  happen  to  coalesce  at  the  location  of 
pressure  measurement  generating  a  high  pressure  region.  This  can  be  seen  from  Figure  5.25 
d  where  the  high  pressure  region  is  concentrated  approximately  at  the  location  of  QSDD 
ABl  whereas  from  Figure  5.25-e  it  can  be  seen  that  the  high  pressure  zone  is  at  the  corners. 
This  phenomena  keeps  on  repeating  itself  resulting  in  the  pressure  peaks  as  can  be  seen  in 
the  pressure  signal  data.  With  the  passage  of  time,  the  strength  of  the  blast  diminishes 
resulting  in  increased  time  difference  between  two  pressure  peaks  and  reduced  pressure  peak 
values. 

Figure  5.26  shows  instantaneous  flow  flelds  of  pressure  for  Room  3  conflguration  from 
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Figure  5.23:  Room  3  with  hall  configuration:  Instantaneous  soot  mass  fraction  at  15.5  and 
34.2  milliseconds 
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Figure  5.24:  Comparison  between  experimental  and  simulation  pressure  history  for  QSDD 
ABl  sensor 
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another  view  (side  view).  The  same  conclusion  can  be  drawn  from  this  figure  as  well. 

Figure  5.27  shows  the  soot  evolntoon  at  this  time  shots.  Untill  2.74  milliseconds  the 
identity  of  the  fireball  can  be  distinguished,  However,  after  that  mixing  dominates  and 
identity  of  the  fireball  is  lost.  Interactions  of  pressure  waves  refiecting  from  the  wall  causes 
disturbances  to  the  fiow  field  and  the  RMI  are  further  enhanced  leading  to  a  more  complex 
mixing  region.  Further  analysis  is  still  needed  to  fully  understand  all  these  features. 

5.4.6  Effect  of  Charge  Shape 

Figure  5.28  shows  effect  of  charge  shape  (cylindrical  vs.  spherical)  on  the  pressure  pulse.  For 
both  the  cases,  mass  and  volume  of  the  charge  has  been  conserved.  The  equivalent  radius  of 
the  spherical  charge  comes  out  to  be  4.404  which  has  same  volume  as  the  cylindrical  charge. 

The  arrival  of  first  pressure  peak  is  off  by  2  milliseconds  for  the  spherical  charge  compared 
to  cylindrical  charge.  Also,  peak  pressure  is  40  psi  less  compared  to  cylindrical  charge.  The 
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Figure  5.25:  Instantaneous  pressure  field  snapshot  for  Room  3  configuration  as  viewed 
from  the  top  at  0,  1,  1.87,  2.74,  3.77  and  4.87  milliseconds 


42 


Pressure  (Pa)  Pressure  (Pa)  Pressure  (Pa) 


Figure  5.26:  Instantaneous  pressure  field  snapshot  for  Room  3  configuration  as  viewed 
from  te  side  at  0,  1,  1.87,  2.74,  3.77  and  4.87  milliseconds 
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Figure  5.27:  Instantaneous  mass  fraction  of  soot  field  for  Room  3  configuration  as  viewed 
from  the  top  and  side  at  1,  1.87,  2.74,  3.77  and  4.87  milliseconds 
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Figure  5.28:  Comparison  between  experimental  and  simulation  pressure  history  for  QSDD 
ABl  sensor  for  C4  cylindrical  vs.  spherical  charge 


characteristics  of  the  initial  pressure  pulses  (till  50  milliseconds)  seem  to  be  quite  different, 
however,  as  time  progresses  the  signals  become  similar.  Figure  5.29  shows  instantaneous 
flow  field  of  pressure  for  the  spherical  charge  case  for  two  different  views  side  and  top 
view.  The  pictures  show  multiple  pressure  wave  coalescence  phenomena  which  is  happening 
intermittently  inside  the  room. 

Figure  5.30  shows  instantaneous  flow  field  of  soot  mass  fraction  for  spherical  charge.  The 
overall  flow  features  are  similar  to  a  cylindrical  charge.  More  analysis  is  stull  underway. 

5.4.7  Effect  of  Charge  Material 

Figure  5.31  shows  effect  of  charge  material  on  the  pressure  pulse.  C4  is  a  stronger  charge 
compared  to  Nitro- methane.  The  peak  pressure  is  110  psi  less  for  NM  charge.  The  time  of 
arrival  of  shocks  is  delayed  for  NM  charge  because  of  decrease  in  shock  velocity. 

Figure  5.32  shows  instantaneous  flow  fields  of  pressure  for  two  different  sectional  views. 
Similar  flow  features  as  compared  to  C4  charge  are  observed  with  Nitro-methane  charge  also. 

Figure  5.33  shows  the  evolution  of  CO  mass  fraction  as  time  progresses.  The  mushroom 
like  structures  caused  by  RM  instability  is  present  in  the  flow  field.  Visibly  these  structures 
are  different  for  C4  and  Nitro-methane  charge. 

5.4.8  Effect  of  After-Burning 

Figure  5.34  shows  effect  of  after  burning  on  the  pressure  pulse.  The  computed  results  for 
reactive  case  match  experimental  data  much  better  than  the  non-reactive  case.  The  peak 
shock  values  and  shock  arrival  time  are  in  better  accordance  with  that  of  the  experiments. 
This  clearly  shows  after-burning  has  huge  impact  on  the  flow  dynamics  of  a  confined  spaced 
environment. 

The  result  shows  higher  value  of  pressure  after  60  milliseconds  compared  to  experiments 
which  might  be  due  to  the  fact  that  this  domain  is  limited  to  the  single  room  configuration. 
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Figure  5.29:  Instantaneous  pressure  field  for  Room  3  configuration  as  viewed  from  the  top 
and  side  at  1,  2.92  and  4.67  milliseconds 
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Figure  5.30:  Instantaneous  mass  fraction  of  soot  field  for  Room  3  configuration  as  viewed 
from  the  top  and  side  at  1,  1.63,  2.92  and  4.67  milliseconds 
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Figure  5.31:  Comparison  between  experimental  and  simulation  pressure  history  for  QSDD 
ABl  sensor  for  C4  vs.  NM  cylindrical  charges 
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Figure  5.32:  Instantaneous  pressure  field  for  Room  3  configuration  as  viewed  from  the  top 
and  side  at  1.2,  3.0  and  5.26  milliseconds 


Figure  5.33:  Instantaneous  CO  mass  fraction  for  Room  3  configuration  as  viewed  from  the 
top  and  side  at  1.2,  3.0  and  5.26  milliseconds 
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Figure  5.34:  Comparison  between  experimental  and  simulation  pressure  history  for  QSDD 
ABl  sensor  for  C4  reactive  and  non- reactive  simulations 


A  full  domain  simulation  is  planned  with  detailed  chemistry  in  future. 

5.4.9  Eulerian-Lagrangian  Results 

As  noted  earlier,  all  the  simulations  have  been  performed  assuming  soot  in  the  gaseous 
phase,  which  is  obviously  erroneous.  In  an  ideal  world,  the  soot  combustion  should  be 
modeled  assuming  it  to  be  present  in  the  solid  phase.  LESLIE  code  uses  E-L  methodology 
to  track  the  particles.  Soot  particles  are  typically  of  3-3.5  micron  in  diameter.  Density  of 
the  soot  is  2200  kg/m3.  A  simple  calculation  will  show  that  to  achieve  soot  mass  fraction 
present  in  te  C4  charge  for  a  20  lb.  charge,  number  of  soot  particles  to  be  tracked  is  2  X  10^^ 
(yes  a  billion  is  just  10®  )!  To  track  these  many  particles  is  beyond  the  hardware  capability. 
Parcel  method  where  every  parcel  contains  identical  number  of  particles  with  same  properties 
might  be  useful  in  this  case.  In  this  study,  an  assumption  has  been  made  as  to  after  the 
detonation  soot  particles  tend  to  agglomerate  and  the  diameter  of  the  particles  increase.  A 
total  of  600000  parcels  are  used  with  100  particles  per  parcel.  With  these  parameters  the 
mass  of  soot  that  can  be  in  dispersed  phase  has  been  calculated  and  remaining  amount  of 
soot  is  kept  in  the  gaseous  phase  as  previously.  This  is  the  dilute  phase  condition  in  E-L 
scenario.  A  preliminary  run  has  been  setup  with  the  above  mentioned  parameters.  A  detailed 
parametric  study  is  planned  in  future.  Figure  5.35  shows  the  particle  evolution  snapshot 
with  time.  Initial  distribution  of  the  particles  is  uniform,  however  as  time  progresses  the 
distribution  becomes  highly  anisotropic.  Effort  is  going  on  to  calculate  the  particle  source 
terms  information  to  quantify  this  anisotropic  behavior. 

Figure  5.36  shows  3D  view  of  the  particle  evolution  snapshot  with  time.  The  particle 
field  shows  the  jetting  pattern  seen  in  unconfined  studies.  Further  studies  are  underway  to 
determine  how  the  post  shock  flow  field  instabilities  such  as  Richtmyer- Meshkov  instability 
and  the  mixing  layer  growth  is  making  these  wavy  and  streaky  patterns. 
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Figure  5.35:  Instantaneous  Soot  particle  trace  colored  by  velocity  for  Room  3  configuration 
as  viewed  from  the  top  at  0.15,  0.40,  0.65,  0.95,  1.25,  1.5,  1.85,  2.15,  2.75  and  3.0  milliseconds 


Figure  5.36:  Instantaneous  Soot  particle  trace  colored  by  velocity  for  Room  3  configuration 
at  0.15,  0.95  and  3.0  milliseconds 
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5.4.10  Spore  Aerosol  Modeling 

In  the  current  study,  droplets  of  a  spore-laden  aqueous  solution  are  introduced  into  the 
domain  of  interest,  i.e.,  re-shock  zone  of  shock  tube  or  in  the  ambiance  of  a  detonated 
explosive  charge.  The  spores  are  considered  to  be  of  the  Bacillus  species.  Although  most 
spores  of  this  species  are  elliptical,  here,  spores  are  assumed  to  be  spherical  with  radius, 
0.4  micron  [14].  The  spore  aerosol  considered  has  a  particular  droplet  size  distribution, 
concentration  of  spores  in  the  initial  solution  and  the  concentration  of  spores  in  the  domain. 
These  parameters  are  set  based  on  past  experimental  studies  [15]  so  that  the  current  results 
can  be  compared  with  the  results  available  in  the  literature. 

When  the  spore  aerosol  interacts  with  hot  gases  in  either  post-shock  or  post-detonation 
flow,  the  water  enveloping  the  spores  evaporates  and  exposes  the  spores  to  the  heat.  Based 
on  the  aerosol  droplet  radius  there  could  be  multiple  spores  in  a  given  droplet.  These  spores 
after  the  evaporation  of  water  can  stay  clustered  or  disperse.  As  the  temperature  of  the 
spores  increases,  based  on  the  quantity  of  heat  received  by  each  spore,  the  spore  kill  can 
occur  due  to  heating  or  mechanical  rupture.  When  spore-laden  aerosol  is  nebulized,  the 
droplets  of  aerosol  are  distributed  in  the  domain  of  interest  with  each  droplet  having  an 
initial  radius  r®.  In  general,  r®  can  be  specified  based  on  a  distribution  function.  However, 
as  the  exact  distribution  of  the  droplet  size  is  not  known,  r®  (microns),  based  on  a  Gaussian 
distribution,  is  specified  as  =  min{ri,  /r^),  where  rj  is  Gaussian  random  variable  with  mean 
Hrj  and  standard  deviation  As  the  values  of  and  an  are  not  available,  a  range  of 
values  are  used  initially  for  the  shock  tube  simulations  and  the  values  which  provide  good 
agreement  with  experimental  results  are  then  used  in  cases  with  explosive  charges.  Also,  the 
initial  atomization  is  assumed  to  reduce  all  droplets  to  size  less  than  Further,  number  of 
spores  per  droplet  is  set  based  on  the  concentration  of  the  spores  in  the  initial  aqueous 
solution.  In  the  current  study,  for  initial  spore  concentration  of  10^0  spores/ml,  varies 
from  1  to  5.  Also,  the  concentration  of  the  spores  in  the  domain,  fin  is  kept  at  1000  spores/cc. 
Due  to  the  heat  transfer  to  the  droplets  in  the  post-shock  region,  the  water  encapsulating 
the  spores  evaporates  and  this  rate  of  mass  transfer  is  given  as  [16] 


m  =  2TrpDrpSh(l  +  Bm)  (5.4) 

Here  D  is  the  diffusivity  of  the  gas.  The  expressions  for  Nusselt  number  (Nu),  Sherwood 
number  (Sh)  and  the  Spalding  mass  transfer  number  {Bm)  for  droplets  are  available  elsewhere 
[16].  After  the  water  evaporates,  no  mass  transfer  is  considered  from  individual  spores  or 
spore  clusters.  Hence,  these  expressions  for  rh  and  Nu  are  used  until  the  radius  of  the  droplet 
reduces  to  the  effective  radius  of  the  spore  cluster  (or  radius  of  the  spore)  present  inside  the 
droplet,  i.e.  .  Here,  the  effective  radius  of  the  spore  cluster  is  determined  as  : 
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Here  fp  is  the  packing  fraction  and  is  taken  to  be  1  when  =  1  and  0.74  otherwise, 
i.e.,  close  packing  assumption.  The  expression  for  Nu  used  when  is  provided  elsewhere  [19]. 
When  the  spores  are  exposed  to  HTGE,  the  temperature  of  the  spores  increase  to  a  critical 
value  {Tc).  Past  calculations  of  thermo- structural  response  of  individual  spores  show  that 
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Figure  5.37:  Initial  Spore  Cloud  location  (marked  by  black  dot)  for  Room  3  configuration, 
Case  A:  Spore  cloud  located  on  the  top  of  the  charge  at  an  angle  45  degree,  Case  B:  Spore 
cloud  located  at  an  angle  45  degree  from  the  charge  in  the  horizontal  plane 


this  critical  temperature  should  result  in  spore  membrane  rupture  and/or  heating  of  spore 
core  leading  to  spore  kill  [17].  Experimental  studies  suggest  that  the  loss  of  spore  viability 
and  structural  damage  occurs  at  gas  temperatures  of  about  750K  and  above  [18].  However, 
for  the  cases  considered  here,  the  exact  quantity  of  heat  needed  to  kill  a  spore  or  to  reach 
Tc  is  not  available.  Hence,  Tc  is  assumed  as  a  variable  parameter  and  a  range  of  values 
are  used  in  each  case  to  obtain  the  percentage  of  spores  killed.  In  cases  presented  here, 
any  spore  whose  temperature  exceeds  Tc  is  assumed  to  be  neutralized.  In  this  study,  spore 
clouds  are  put  at  different  locations  of  the  room  with  different  orientation  from  the  charge 
and  thereafter  spore  dispersion  has  been  calculated.  Figure  5.37  shows  the  initial  spore  cloud 
location  for  these  various  cases.  Two  different  cases  have  been  chosen  for  the  preliminary 
study  Case:  A  has  the  spore  cloud  located  on  the  top  of  the  charge  at  45  degree  angle. 
Case:  B  has  the  spore  cloud  near  to  the  wall  but  at  a  45  degree  angle. 

Figure  5.38  shows  the  spore  dispersion  for  Case  A  with  time. 

For  initial  5  milliseconds  the  spore  cloud  is  not  affected  by  the  detonation  as  shock  wave 
does  not  hit  it.  Once  the  shock  passes  through  the  cloud,  the  cloud  starts  to  disperse  and  at 
the  same  time  gains  temperature  after  interacting  with  the  shock.  Dispersion  of  the  spore 
cloud  increases  with  time  leading  to  a  variation  in  the  spore  temperature  as  the  spores  mix 
in  the  post  detonation  environment.  Further  studies  of  spore  transport  and  the  entrainment 
into  the  turbulent  mixing  layers  are  needed  to  understand  the  dynamics  but  such  studies 
requires  extensive  data  processing  and  these  tools  are  being  created  for  this  purpose. 

Figure  5.39  shows  the  spore  survival  rate  as  a  function  of  time  for  the  two  different  cases 
mentioned  above.  Clearly  being  close  to  the  charge  helps  in  increasing  post  detonation  spore 
kill  but  in  general  there  are  other  issues  to  consider.  Preliminary  analysis  also  shows  that 
this  is  also  related  to  the  entrainment  in  the  mixing  layers  behind  the  detonation  front. 
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Figure  5.38:  Evolution  of  Spore  Cloud  (marked  by  black  dots)  for  Room  3  configuration, 
for  Case  A  at  6,  16.5,  20.6  and  24.5  milliseconds 
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Figure  5.39:  Percentage  of  spores  left  intact  as  a  function  of  time  for  Case  A  and  Case  B 
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5.5  Detonation  in  a  Two  Room  Configuration 

5.5.1  Simulation  Setup 

To  investigate  the  interaction  of  spores  with  hot  detonation  product  gases  in  a  confined  envi¬ 
ronment,  a  two  room  configuration,  shown  in  Fig.  5.40  is  used.  The  two  room  configuration 
comprises  of  a  inner  room  enclosed  by  an  outer  room.  The  rooms  are  connected  by  a  vent 
hole  of  diameter  0.1  m.  The  length  and  the  width  of  the  inner  room  are  3  m  each,  and  the 
height  is  1  m.  The  outer  room  is  9  m  in  width  and  length,  and  5  m  in  height.  The  rooms  are 
separated  by  a  wall  of  0.5  m  thickness.  The  vent  hole  is  located  at  the  center  of  roof  of  the 
inner  room.  All  the  walls  are  assumed  to  be  adiabatic  free-slip  walls.  A  cylindrical  explosive 
is  placed  on  fioor  of  the  inner  room  directly  below  the  vent  hole.  The  explosive  charge  is 
22.86  cm  in  height  and  2.286  cm  in  radius.  Small  perturbations  are  introduced  on  the  initial 
charge  surface  to  mimic  the  natural  surface  imperfections  and  induce  hydrodynamic  insta¬ 
bilities.  The  two-room  configuration  is  discretized  using  44  million  computational  cells.  The 
minimum  size  of  the  computational  cell  is  set  to  0.6  mm  and  the  cells  are  distributed  in  the 
domain  to  resolve  all  the  relevant  flow  features. 

The  explosive  material  composed  of  HMX  with  inert  metal  particles.  The  number  of 
metal  particles  in  the  charge  is  0.5  million,  and  the  radius  of  each  particle  is  20  micrometers. 
Therefore,  the  volume  fraction  of  charge  is  about  4.46  x  10“^.  In  addition,  the  0.1  million 
spore  particles  are  placed  at  the  corner  of  the  inner  room.  They  are  also  treated  as  inert 
particles.  The  radius  of  each  spore  is  30  micrometers. 

Unlike  the  procedure  developed  earlier  [73,  74],  the  initial  blast  wave  is  not  modeled 
based  on  the  one-dimensional  detonation  profile  but  is  generated  by  the  full  three-dimensional 
detonation  propagation  through  the  condensed  phase  charge.  To  achieve  this,  a  high  pressure 
zone  is  placed  at  the  center  of  the  explosive  charge  and  a  detonation  is  allowed  to  propagate 
in  axial  direction  of  the  charge.  The  detonation  is  is  expected  to  reach  the  top  and  the 
bottom  of  the  charge.  Thus,  the  simulation  accounts  for  a  wide  range  of  length  and  time 
scales  from  the  detonation  (small  scale)  to  the  blast  wave  propagation  (large  scale). 

An  Eulerian-Lagrangian  method  is  used  to  solve  the  governing  equations  [73,  75].  A 
combination  of  Mie  Gruneisen,  Jones-Wilkins-Lee  (JWL),  and  thermally  perfect  gas  equa¬ 
tions  of  state  (EOS)  are  employed;  the  Mie  Gruneisen  EOS  is  applied  in  the  regions  with 
condensed  phase  charge,  while  Jones-Milkins-Lee  EOS  and  thermally  perfect  EOS  are  used 
for  gas  phase  detonation  products  and  ambient  air. 

5.5.2  Results  and  Discussion 

The  chronology  of  evolution  of  the  blast  wave  and  the  post-detonation  mixing  zone  are  shown 
in  Fig.  5.41.  The  detonation  of  the  initial  charge,  at  t  =  0.0,  generates  a  high  pressure  and 
a  high  temperature  zone  at  the  center  of  the  charge.  The  region  of  high  temperature  spreads 
in  the  axial  direction  behind  the  detonation.  The  hot  gaseous  products  formed  vent  out  in 
the  radial  direction  out  of  the  charge.  Also,  the  particles  in  the  charge  are  ejected  but  remain 
within  the  hot  gases.  At  this  stage,  the  particles  acquire  momentum  and  energy  from  the 
detonation  product  gases  but  due  to  inertia  are  slower  than  the  detonation  products. 

At  about  20/is,  the  condensed  phase  explosive  is  completely  consumed,  resulting  in  the 
formation  of  the  primary  blast  wave.  The  blast  wave  propagating  towards  the  ground,  reflects 
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(a)  Full  view  of  the  two  rooms 


(b)  Enlarged  view  of  the  inner  room 

Figure  5.40:  Geometry  used  for  this  simulation.  Black  cylinder:  explosives,  blue  sphere: 
Spore  particle,  green  cylinder:  Vent  hole  connecting  inner  room  and  outer  room 
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from  the  floor,  and  produces  a  wrinkled  flreball  near  the  floor  (Fig. 5. 41  (c)).  In  addition,  the 
particles  penetrate  the  contact  surface  between  the  detonation  products  and  the  ambient  air. 
As  the  particles  overtake  the  contact,  the  particles  induce  additional  perturbations  i.e.,  in 
addition  to  that  introduced  by  the  imperfections  on  the  initial  charge  surface.  They  are  not 
signiflcant  in  this  phase,  but  affect  the  hydrodynamic  instabilities,  as  is  observed,  at  later 
stages. 

At  t  =  0.3  ms,  the  blast  wave  eventually  reaches  the  ceiling  of  the  inner  room  and 
enters  the  outer  room  through  the  vent  hole  (Fig. 5.41  (e)).  The  particles  entrained  into  the 
mixing  zone  are  also  ejected  into  the  outer  room.  The  blast  wave  and  the  gaseous  detonation 
products  exhausted  from  the  vent  hole  form  a  mushroom  shaped  structure(Fig.5.41  (f)).  The 
constriction  by  the  vent  hole  and  the  strength  of  the  initial  blast  do  not  permit  the  flreball 
to  expand  into  the  outer  room  beyond  the  conflnes  of  the  vent  hole.  In  inner  room,  the 
temperature  increases  from  0.3  to  0.9  ms  due  to  the  compression  wave  generated  to  recover 
the  over-expanded  region  at  center  of  the  inner  room.  Also,  during  this  phase,  the  small 
perturbations  that  is  introduced  initially  or  added  by  penetrations  of  particles  develop  and 
add  signiflcant  asymmetric  to  the  shape  of  the  flreball.  They  are  especially  signiflcant  at  the 
upper  portion  of  the  blast  wave  and  the  flreball. 

Around  0.9  ms,  the  primary  blast  wave  hits  and  reflects  from  side  wall/ceiling.  The 
reflected  blast  waves  interacts  with  the  detonation  products  and  form  complex  hydrodynamic 
structures.  In  particular,  the  hydrodynamic  instabilities  due  to  RMI  and  RTI  occur  at  the 
location  at  which  the  reflected  shock  wave  encounters  the  surface  of  detonation  products. 
When  the  blast  wave  contacts  to  the  spore  particles,  spore  particle  cloud  is  deformed  along 
the  contact  surface  of  the  shock.  The  spores  are  driven  towards  the  wall,  and  pushed  onto 
the  corner  of  the  inner  room  (Fig. 5.41  (h))  with  relatively  higher  gaseous  temperature. 


5.6  Multi-Blast  Wave  Interactions 

Interaction  of  multiple  blast  waves  can  be  used  to  direct  energy  toward  a  target  while  si¬ 
multaneously  reducing  collateral  damage  away  from  the  target  area.  In  [3] ,  authors  simulate 
multiple  point  source  explosives  and  the  resulting  shock  interaction  and  coalescence  behavior 
were  explored.  Different  munitions  were  placed  concentrically  around  the  target.  Different 
patterns  are  found  and  particularly,  a  strong  linear  behaviour  is  observed. 

5.6.1  Two  Blast  Wave  Sources 

To  study  the  interaction  of  two  blast  waves,  we  mimic  the  set  up  from  [3].  In  this  study, 
researchers  used  a  parrallel  flnite  difference  solver  on  overlapping  grids  to  solve  the  Euler 
equations  with  shock  capturing  second  order  Godunov  Scheme.  They  use  adaptative  mesh 
refinement  method  to  reduce  the  CPU  cost  of  the  simulations.  Taylor  similarity  laws  [76] 
for  pressure,  density  and  velocity  are  used  to  initialize  the  blast  spot.  Different  setups 
are  presented  to  show  interaction  of  blast  waves  in  a  2D  domain.  In  particular,  the  study 
focused  on  the  interaction  af  blast  waves  placed  around  a  target  point,  when  all  the  munitions 
detonate  initially  at  the  same  time.  In  the  present  study,  the  similar  simulations  are  going 
to  be  compared  to  these  reference  cases. 
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(a)  t  =  Qjis 


(b)  t  =  10/is 


0.000050  s  0.000050  s 


(c)  t  =  50/is 


(d)  t  =  0.1ms 
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(e)  t  =  0.3m5 


(f)  t  =  0.6ms 


(g)  t  =  0.9ms 


(h)  t  =  1.3ms 


Figure  5.41:  Formation  and  Evolution  of  blast  wave.  Left:  Particles  dispersion  with  iso¬ 
surface  of  detonation  products.  Right:  Temperature  profiles(min:  300  K,  max:  2100  K.) 
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Figure  5.42:  Taylor  profiles  for  pressure,  density  and  velocity  :  our  current  simulation 
initialization  is  similar  with  the  reference  one 


The  Taylor  similarity  laws  are  used  to  initialize  the  profiles  of  pressure,  density  and 
velocity  inside  the  blast  wave  initially.  A  comparison  of  the  profiles  used  in  the  reference 
simulations  by  [3]  and  the  current  simulations  are  shown  in  Fig.  5.42.  While  the  velocity 
behaves  quasi  linearly  inside  the  blob,  density  and  pressure  have  an  exponential  trend,  with 
maximnm  valne  at  the  external  edge  of  the  blast  wave.  Note  that  our  current  simulations 
use  a  similar  initialization  compared  to  the  reference  ones. 

The  evolution  of  two  blast  waves  is  necessarily  the  same  for  both  waves  until  both  pressure 
front  interact.  If  the  shock  is  weak,  the  interaction  is  pretty  straightforward,  and  can  be 
schematised  with  Fig.  5.43  a).  Indeed,  in  this  case,  both  waves  behave  like  they  do  not 
see  each  other.  Note  that  this  no  interaction  case  has  been  simnlated  in  the  present  work, 
but  the  results  are  not  shown.  If  the  pressure  front  is  strong  enough,  non  linear  interaction 
occurs  (Fig.  5.43  b)).  At  this  point,  the  spherical  shape  of  the  blast  waves  is  lost.  This 
result  has  been  reproduced  in  the  current  simulation  in  Fig.  5.43  c).  In  this  case,  the  oval 
interaction  zone  is  not  present,  and  a  similar  diamond  shape  appeared. 

5.6.2  Blast  wave  inside  a  box 

The  ontfiow  conditions  that  were  nsed  to  avoid  reflections  at  the  bonndaries  are  now  removed, 
and  no  slip  walls  are  used.  The  goal  is  to  see  the  reflection  of  the  blast  wave  at  the  boundary. 
The  initialisation  is  exactly  the  same  in  the  case:  pressure,  density  and  velocity  profiles  are 
specified  inside  a  1.5  mm  blob  inside  the  2D  domain.  Pressure  fields  at  different  times  are 
presented  in  Fig. 5. 45.  Note  that  as  the  pressure  amplitude  is  decreasing  with  time,  each 
pressure  field  has  been  rescaled  to  see  the  pressure  front. 
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(b)  Simulation  from  [3]  (c)  Current  simulation 


Figure  5.43:  Qualitative  comparison  of  two  blast  waves  interaction  :  Schematic  of  2  blasts 
wave  interaction,  comparison  between  the  simulation  from  [3]  and  our  simulation  :  our  work 
is  able  to  capture  a  non  linear  behaviour 


(a)  From  [3] 


(b)  Actual  profiles 


Figure  5.44:  Pressure  profile  from  the  blast  wave  origins 
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(g)  At  t=  53.4  fis  (h)  At  t=  80.3  fis  (i)  At  t=  133.3  /is 

Figure  5.45:  Pressure  fields  with  time  :  the  blast  wave  evolves  with  a  spherical  shape  until 
it  reaches  a  wall.  At  this  point,  reflection  occnrs,  and  interaction  happens 
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Figure  5.46:  Speed  up  for  the  LESLIE  Solver 


From  Fig. 5. 45  a)  to  c),  the  evolution  of  the  blast  wave  is  very  classical  :  the  blast  wave 
evolves  with  a  spherical  shape.  Then  a  part  of  the  pressure  contour  reaches  the  boundary, 
which  reflects  the  blast  wave  inside  the  domain  in  Fig. 5. 45  d)  and  e).  At  this  point,  the  shape 
of  the  shock  front  is  still  spherical.  Then,  a  non  linear  behaviour  occurs  :  the  spherical  shape 
of  the  blast  wave  disappeared,  and  a  cusp  point  appeared  at  the  center  (Fig. 5. 45  f)  and  g)). 
This  modified  blast  wave  interacts  then  with  a  secondary  blast  wave  (Fig.  5.45  h)).  Finally, 
the  primary  blast  wave  interacts  with  the  other  walls  and  corners  to  create  a  much  more 
complex  shape. 

5.7  Performance  and  Scaling  Analysis 

As  explained  earlier  in  this  report,  the  Lagrangian  simulation  for  dense  clouds  can  be  very 
expensive  in  terms  of  CPU  cost  because  of  the  high  number  of  particles  it  needs  to  solves. 
That  is  why  the  use  of  an  Eulerian  solver  for  the  dense  phase  coupled  with  the  Lagrangian 
Solver  for  the  dilute  phase  is  crucial  for  a  big  range  of  mass  loading.  For  the  Lagrangian 
Solver,  Fig. 5. 46  shows  the  speed  up  properties  of  the  code  for  different  numerical  schemes 
and  different  number  of  particles.  Performance  drops  when  the  number  of  particles  increases. 
Note  that  Mac  means  MacCormack  scheme. 
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5.8  Multi-scale  Modeling  of  Condensed  Phase  Detona¬ 
tion 

Blast  waves  are  typically  initialized  using  well-known  one-dimensional  detonation  profiles 
that  have  been  validated  against  experimental  observations.  While  this  assumption  is  valid 
for  cases  of  one-dimensional  nature,  for  example  the  explosion  of  a  spherical  homogeneous 
charge,  it  does  not  account  for  properties  pertaining  to  the  condensed  phase  energetic  ma¬ 
terial.  The  shape  of  the  explosive,  as  well  as  its  structural  and  chemical  composition,  can 
have  significant  implications  on  the  characteristics  of  the  formed  blast  wave,  including  its 
speed,  dispersion,  and  level  of  instability.  For  example,  a  randomly  packed  energetic  mate¬ 
rial  develop  turbulent  breakdown  early  in  the  process  due  to  the  seemingly  perturbed  energy 
depositions.  Variations  at  the  stage  of  the  transition  to  detonation  can  develop  and  grow  to 
produce  outcomes  that  would  otherwise  be  unexpected.  Recent  advancements  in  the  field 
enable  such  simulations  to  be  initialized  with  the  unburnt  condensed  phase  material  to  bet¬ 
ter  predict  the  blast  wave[77,  78].  The  procedure  introduces  the  advantage  of  quantifying 
uncertainties  in  terms  of  the  explosive  content,  binder  material,  and  granularity.  By  doing 
so,  the  approach  better  predicts  blast  wave  characteristics,  thus  isolating  the  uncertainty 
quantification  of  the  numerical  models  while  taking  into  account  the  stochastic  nature  of  the 
physical  problem. 

In  summary,  it  is  necessary  to  address  this  multi-scale  problem  starting  with  the  unburnt 
condensed  phase  explosive  for  the  following  3  main  reasons: 

1.  The  simulation  better  predicts  the  blast  wave  characteristics,  such  as  speed,  shape, 
and  dispersion,  based  on  the  explosive’s  properties. 

2.  The  perturbations  at  the  micro-scale  develop  early  turbulent  breakdown  and  mixing 
that  may  have  significant  implications  on  the  outcome. 

3.  The  results  are  based  on  the  stochastic  nature  of  the  condensed  phase  enabling  a  better 
quantification  of  the  uncertainties  arising  from  the  numerical  modeling. 

5.8.1  Formulation  of  the  Equation  of  State 

The  Mie-Gruneisen  equation  of  state  (EOS)  is  generally  accepted  to  perform  well  with  solid 
energetic  materials  [79,  80].  Its  validity  extends  beyond  the  detonation  wave  front  and 
reaction  zone  to  the  point  where  the  expanding  products  reach  a  lower  density  at  which 
the  Mie-Gruneisen  EOS  starts  to  deteriorate  [79].  At  this  point,  another  equation  of  state, 
which  is  introduced  in  a  later  section,  is  needed  for  the  proper  resolution  of  the  evolution  of 
the  gaseous  products.  The  Mie-Gruneisen  EOS  is  expressed  as: 

e{P,i^)  =  ^[P  -  f  (z^)]  +  eo  (5.6) 
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The  specific  volume  and  the  pressure  along  the  Hugoniot  are  herein  represented  by  u  and 
P-H,  and  s  is  defined  as  the  slope  in  the  relationship  between  the  shock  speed  and  the  particle 
velocity,  where  Ug  =  cq  +  sUp. 

To  allow  for  the  formation  of  a  blast  wave,  the  detonation  front  propagates  the  entire 
condensed  domain  and  gets  released  into  a  gaseous  surrounding  atmosphere,  Air.  One  of 
the  challenges  facing  the  complete  simulation  of  such  problems  is  the  proper  handling  of  the 
thermodynamic  properties.  On  one  hand,  the  Mie-Gruneisen  equation  of  state  is  successfully 
used  to  handle  the  high  density  phase  of  the  domain.  On  the  other  hand,  it  rapidly  deterio¬ 
rates  as  the  density  of  the  expanding  detonation  products  falls  below  the  reference  density 
used,  Pq.  For  the  surrounding  air  at  atmospheric  conditions,  it  can  be  treated  using  the  ideal 
gas  equations.  As  for  the  stage  in  between,  the  Jones-Wilkins-Lee  (JWL)  equation  of  state 
seems  to  be  adequate  as  it  has  been  developed  for  such  purposes.  The  JWL  equation  of  state 
has  been  used  in  numerous  studies  to  calculate  the  thermodynamic  properties  of  detonation 
product.  It  is  expressed  as: 
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The  complete  simulation  of  the  shock  to  detonation  transition  until  after  the  formation  of 
a  blast  wave  requires  the  use  of  all  three  equations  of  state.  The  combination  of  the  equations 
of  state  is  implemented  in  terms  of  the  density  of  the  products,  where  density  limits  are 
specified  to  denote  the  regions  of  transition.  For  expanding  products  with  densities  higher 
than  the  upper  limit  of  the  condensed  phase  pcu^  the  Mie-Gruneisen  equation  of  state  is 
used,  whereas  JWL  is  used  when  it  drops  below  the  lower  limit  pcL-  A  smoothing  function 
is  used  for  the  range  in  between.  A  similar  procedure  is  used  for  the  transition  between  JWL 
and  ideal  gas  equations  using  upper  and  lower  density  limits  pou  and  pcL-  The  resulting 
expression  takes  the  following  form: 


/ 


X  =  < 


(a^MG  —  3;jwl) 
a^jwL 

(3:^  JWL  —  3:ig) 


p  —  PCL 
Pcu  —  PCL 

P  —  Pgl 

pGU  —  pGL 


P  >  Pcu 
Pgl  <  P  <  Pgu 
Pgu  <  P  <  Pgl 
Pgl  <  P  <  Pgu 


(5.9) 


I,  3^10  p  <  Pgu 

where  x  represents  thermodynamic  variables.  Figure  5.47  compares  P  —  u  plots  for  the  three 
equations  of  state  alongside  the  combined  equation  for  HMX  at  5000  K.  The  density  limits  are 
chosen  in  such  a  way  to  smoothen  the  transition.  For  this  case,  pcu,  Pcl,  Pgu,  and  pgl  are 
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Figure  5.47:  P  —  p  plot  for  the  combined  equation  of  state  with  respect  to  the  Mie- 
Griineisen,  JWL  and  ideal  gas  equations  of  state. 

2500  kg/m^,  1000  kg/m^,  1000  kg/m^,  and  500  kg/m^,  respectively.  The  transition  from  the 
condensed  phase  to  the  expanded  blast  wave  products  is  very  rapid,  and  therefore,  variations 
in  the  density  limits  have  negligible  effects  on  the  resulting  blast  wave  characteristics.  Figure 
5.47  shows  a  smoothened  transition  from  the  MG  to  IG  through  the  JWL  equation  of  state. 

5.8.2  Simulations  Setup 

The  purpose  of  this  study  is  to  investigate  the  effects  of  the  micro-structure  of  the  condensed 
phase  explosive  on  the  resulting  blast  wave.  The  simulations  consider  modifications  in  the 
condensed  phase  that  are  typically  ignored  in  blast  wave  studies.  Three  cases  are  addressed 
using  charges  in  the  shape  of  a  flat  plate.  The  shape  introduces  two-dimensionality  effects 
that  can  change  the  dynamics  of  the  blast  wave,  as  seen  in  experiments  [81].  The  first 
case  consists  of  a  homogeneous  HMX  charge  being  initiated  with  a  strong  cylindrical  shock 
through  its  middle  section,  as  shown  in  Fig.  5.48(a).  The  shock  represents  an  initiation 
similar  to  that  produced  by  a  detonator.  The  yellow  color  denotes  the  energetic  material 
and  the  blue  part  is  air.  The  red  half-circle  shows  the  location  of  the  initiating  shock.  In 
Fig.  5.48(b),  the  second  case  uses  the  same  setup  but  using  a  Polymer-Bonded  Explosive 
(PBX)  with  randomly  packed  HMX  crystals  in  an  Estane  binder.  Being  more  realistic,  the 
case  includes  the  natural  perturbations  in  the  energetic  material  to  study  their  effects  on  the 
blast  wave  formation.  The  last  case  also  makes  use  of  the  same  PBX,  however,  the  initiating 
shock  is  located  on  the  top  portion  of  the  charge  at  a  part  in  contact  with  air.  Moreover, 
the  shock  strength  is  tremendously  reduced,  and  uses  a  hot  spot  to  aid  in  the  ignition,  and 
it  is  referred  to  as  the  weak-shock  case. 

5.8.3  Results  and  Discussion 

Eigure  5.49  shows  the  velocity  contour  plots  at  3.6  (is.  At  this  point,  detonation  waves  form 
and  propagate  through  the  energetic  material.  When  comparing  the  3  cases,  the  detonation 
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(a)  Homogeneous  HMX  (b)  PBX-Strong  shock  (c)  PBX-Weak  shock 


Figure  5.48:  Diagrams  of  the  initial  setup.  The  HMX  material  is  represented  by  the  yellow 
color  and  Air  is  in  blue.  In  the  cases  with  PBX,  the  HMX  cystals  are  embedded  in  Estane 
as  the  binder. 


(a)  Homogeneous  HMX  (b)  PBX-Strong  shock  (c)  PBX-Weak  shock 


Figure  5.49:  Shape  of  the  detonation  waves  at  3.6  /ts.  The  PBX  cases  show  more  pertur¬ 
bations  in  the  wave  due  to  the  irregularities  in  the  HMX  crystals  packing. 

wave  in  the  homogeneous  case  has  traveled  a  bigger  portion  of  the  charge.  It  has  made  its 
way  out  of  the  material  at  the  middle  top  part  where  the  blast  wave  starts  forming.  The 
overall  shape  of  the  wave  assumes  a  smooth  edge  all  around.  Looking  at  the  PBX  case  in 
Fig.  5.49(b),  the  detonation  front  is  irregular,  and  has  yet  to  travel  the  height  of  the  charge. 
As  for  the  weak-shock  case,  the  blast  wave  starts  forming  as  the  products  eject  into  air  while 
the  detonation  wave  burns  through  the  charge.  When  compared  to  the  other  cases,  the 
amount  of  the  burnt  part  of  the  slab  is  noticeably  less. 

The  velocity  contour  plots  for  the  3  cases  are  shown  in  Fig.  5.50  at  46  fis  at  which  time  the 
material  is  completely  burnt.  The  wave  front  of  the  homogeneous  HMX  has  progressed  more 
than  the  other  cases.  This  is  expected  since  the  binder  not  only  occupies  room  in  the  slab 
thus  reducing  the  total  amount  of  HMX,  but  its  reaction  is  also  endothermic,  which  draws 
energy  from  the  reacted  energetic  material .  Figure  5.50(a)  also  shows  a  smooth  symmetrical 
structure  of  the  propagation.  While  for  the  PBX  cases,  an  irregular  and  asymmetric  behavior 
forms. 

Looking  at  the  progress  of  the  fireball  in  Fig.  5.51,  the  expansion  of  products  fall  behind 
the  blast  wave  front.  The  homogeneous  case  demonstrates  a  symmetric  expansion,  similar 
to  what  is  routinely  used  in  the  simulation  of  blast  waves.  However,  when  examining  the 
PBX  cases,  an  asymmetric  random  structure  emerges.  The  front  of  the  product  expansion 
exhibits  finger- like  instabilities  as  seen  in  experiments.  Its  justification  stems  from  the  initial 
random  distribution  of  crystals  at  the  micro-structure.  It  is  interesting  to  note  that  although 
the  weak-shock  case  is  the  first  to  develop  a  blast  wave  due  to  the  location  of  the  initiating 
shock,  its  progress  is  the  most  retarded.  The  tardiness  of  the  PBX  cases  is  better  justified 
when  examining  the  degree  of  mixing  and  CO2  generation  in  a  later  discussion. 
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Velocity  (m/s]  Magnitude 


Velocity  (m/s)  Magnitude 


Velocity  (m/s)  Magnitude 


(a)  Homogeneous  HMX  (b)  PBX-Strong  shock  (c)  PBX-Weak  shock 

Figure  5.50:  Shape  of  the  blast  waves  at  46  fis.  The  blast  from  the  homogeneous  case  is 
at  a  more  advanced  stage  while  the  PBX  cases  show  more  irregular  shapes. 
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(a)  Homogeneous  HMX  (b)  PBX-Strong  shock 


(c)  PBX-Weak  shock 


Figure  5.51:  Dispersion  of  products  at  46  /is.  The  PBX  cases  show  more  prominent 
instability  structures  due  to  the  irregularities  in  the  initial  packing  at  the  micro-scale. 
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(a)  Homogeneous  HMX  (b)  PBX-Strong  shock  (c)  PBX-Weak  shock 


Figure  5.52:  CO2  mass  fraction  at  46  /rs.  Although  the  homogeneous  case  is  at  a  more 
advanced  stage,  the  higher  generation  of  CO2  in  the  PBX  cases  is  indicative  of  the  higher 
rate  of  mixing. 


Two-dimensional  plots  of  the  CO2  mass  fraction  are  provided  in  Fig.  5.52.  The  generation 
of  carbon  dioxide  is  the  result  of  the  afterburn,  such  that: 

2CO  +  02^  ‘2CO2 

Since  its  generation  is  the  outcome  of  the  interaction  of  the  reaction  products  with  air,  the 
CO2  mass  fraction  can  be  a  good  indicative  of  the  level  of  mixing  and  turbulence  of  the 
expanding  cloud.  The  initial  perturbations  in  the  PBX  cases  stir  the  flow  to  develop  more 
mixing  and  therefore  generate  more  CO2.  The  after-burn  takes  place  at  the  front  of  the 
expanding  products  which  can  be  the  reason  behind  the  name  “fireball” . 

The  degree  of  mixing  is  a  good  measure  of  the  level  of  turbulence  and  mixing  of  the 
reaction  products  with  the  surrounding  air.  For  this  simulation,  the  degree  of  mixing  can  be 
defined  by: 

Degree  of  Mixing  =  (hA.  +  JcoJ  dV 

»^reball  J  ^H20  dV 

Figure  5.53(a)  compares  the  time  evolution  of  the  degree  of  mixing  of  the  3  cases.  In  the 
initial  phases,  it  may  seem  that  the  weak-shock  case  has  the  most  mixing,  however,  this  is 
because  of  the  location  of  the  initiating  shock  giving  the  ejecta  a  head  start  over  the  other 
cases.  While  the  PBX  with  strong  shock  is  the  latest  to  start  mixing  with  air,  the  higher 
speeds  along  with  the  initial  perturbations  develop  stronger  mixing. 

Another  aspect  of  interest  is  a  length-scale  that  compares  the  progress  of  the  blast  wave 
to  that  of  the  expanding  products.  This  length-scale  is  a  measure  of  the  distance  between 
the  blast  wave  front  and  products  front  and  is  labeled  here  as  the  blast  wave  thickness. 
One  way  to  compute  it,  and  to  work  around  the  asymmetric  behavior,  is  to  compute  the 
volume  covered  by  the  blast  wave  and  subtract  the  volume  covered  by  the  reaction  products. 
Plotted  in  Fig.  5.53(b),  the  higher  thickness  signifies  a  faster  progress  of  the  blast  wave  over 
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(a)  Degree  of  Mixing  (b)  Blast  wave  thickness  (c)  Generation  of  C02 


Figure  5.53:  Comparison  of  the  evolution  of  a)  Degree  of  mixing,  b)  blast  wave  thickness, 
and  c)  CO2  generation  with  time  between  the  homogeneous  HMX,  PBX  with  strong  shock, 
and  PBX  with  weak  shock  initiation. 

the  products.  The  comparison  shows  that  the  case  with  PBX  and  strong  shock  initiation 
has  the  largest  thickness  as  the  flow  progresses.  This  is  attributed  to  the  larger  degree  of 
mixing  seen  by  these  products  which  contributes  to  a  slower  expansion.  As  for  the  PBX 
with  weak-shock  initiation,  the  initial  jump  is  the  progress  of  a  shock  rather  than  a  blast 
wave  which  has  yet  to  reach  the  hot  spot  and  ignite.  However,  after  ignition,  the  progress  of 
the  product  expansion  is  closer  to  the  blast  wave  front,  when  compared  to  the  other  cases, 
possibly  due  to  the  lower  speeds  that  this  case  is  undergoing. 

Lastly,  a  comparison  of  the  time  evolution  of  the  CO2  generation  is  plotted  in  Fig.  5.53(c). 
The  second  case  shows  a  significantly  higher  rate  of  carbon  dioxide  generation.  Also  related 
to  the  higher  levels  of  mixing,  the  initial  random  packing,  along  with  the  strong  shock, 
develop  a  turbulent  breakdown  earlier  in  the  simulation.  Although  the  homogeneous  HMX 
case  has  progressed  the  most  at  this  time,  it  has  the  least  amount  of  instabilities  and  therefore 
CO2  generated. 
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CHAPTER  VI 


CONCLUSIONS  AND  METRICS 


6.1  Conclusions 

Under  this  project  a  new  multi-scale  simulation  capability  was  developed  and  demonstrated 
for  varied  applications  of  DTRA’s  interest.  It  is  noted  that  this  effort  required  leveraging 
other  funding  and  man  power  funded  under  Basic  Rsearch  6.1  programs  of  DTRA  and  ONR 
but  the  combined  effort  has  established  a  new  simulation  tool  that  can  be  deployed  for 
many  complex  problems.  Preliminary  validation  of  the  two-phase  solver  and  its  ability  are 
reported  in  this  final  report.  Additionally,  some  of  the  subgrid  models  were  transitioned 
to  our  collaborators  under  this  FRBAA  (CRAFT  Tech)  and  they  have  successfully  applied 
their  code  to  problems  of  specific  interest  to  DTRA.  A  detailed  report  summarizing  their 
effort  has  already  been  submitted  to  DTRA  in  early  2016. 

The  multi-scale  simulation  capability  developed  in  this  effort  is  being  further  applied  to 
other  applications  of  related  interest.  The  predictive  capability  of  this  solver  still  requires 
further  investigation  by  applying  it  to  new  problems  and  such  studies  will  be  considered  in 
the  near  future. 


6.2  Metrics 

6.2.1  Publications 

The  following  publications  were  supported 

1.  Menon,  S.  and  Gottiparthi,  K.  C.,  “Detonation  in  Multi-room  Structure,”  Video  pro¬ 
vided  to  DTRA  for  presentation,  Aug  2015. 

2.  Fedina,  K,  Gotiparthi,  K.  C.,  Fureby,  C,  and  Menon,  S.,  “Combustion  in  Afterburning 
Behind  Explosive  Blasts,”  in  Coarse  Grained  Turbulent  Mixing  (Fernando,  E.,  Ed.), 
Gambridge  University  Press,  2016  (in  press). 

3.  Gottiparthi,  K.  G.,  Schulz,  J.  G.,  and  Menon,  S.,  “On  the  Neutralization  of  Bacterial 
Spores  in  Post  Detonation  Flows,”  Shock  Waves,  Vol.  24,  pp.  455-466,  2014. 

4.  Gottiparthi,  K.  C.,  Schulz,  J.  C.,  and  Menon,  S.,  “Uncertainty  Quantification  of  Bac¬ 
terial  Aerosol  Neutralization  in  Shock  Heated  Gases,”  Shock  Waves,  Vol.  25,  pp.  pp. 
77-90,  2015. 

5.  Menon,  S.,  “Turbulent  Mixing  and  Afternburn  in  Post-Detonation  Flow  with  Dense 
Particle  Glouds.”  Invited  Paper,  19th  Biennial  Gonference  of  the  APS  Topical  Group 
on  Shock  Gompression  of  Gondensed  Matter,  June  2015  (to  appear,  2016). 
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6.2.2  Students  and  Staff 

During  the  course  of  this  project  the  following  people  were  supported  in  Georgia  Tech  dnring 
the  period  Jnne  2011  -  June  2016: 

1.  Post  Doctoral  Fellows:  2  (partial  support  over  2-3  years) 

•  Dr.  Michel  Akiki 

•  Dr.  Gregory  Hannebique 

2.  Ph.D.  Students:  2  (co-shared  with  another  DTRA  project) 

•  Joseph  Schulz,  2015:  A  Study  of  Magnetohydrodynamic  effects  in  turbulent  su¬ 
personic  flows  with  application  to  detonation  and  explosion 

•  Kalyana  Gottiparthi,  2015:  A  Study  of  Dispersion  and  Combustion  of  Particle 
Clouds  in  Post  Detonation  Flows 

3.  M.S.  Students  (without  thesis):  1 

•  Yusuke  Nagaoke,  July  2015 

4.  Undergraduate  Student  Special  Topic:  2 
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1.0  YEAR  1  SUMMARY 


1.1  Project  Major  Goals 

The  primary  goals  for  the  Year  1  work  period  under  this  subeontraet  effort  were  as  follows: 

1  Assemble  a  database  of  BW  agent  neutralization  models. 

2  Develop  a  particle  -  turbulence  interaction  model  for  application  to  the  spore 
neutralization  model. 

1.2  Accomplishments 

During  the  Year  1  work  period,  efforts  focused  on  Tasks  1.2,  1.4  and  1.5.  Under  Task  1.2,  a 
literature  review  of  current  state-of-the-art  biological  warfare  neutralization  models  was 
completed.  This  survey  identified  the  form  of  the  neutralization  model  and  identified  the  model 
correlation  coefficients  for  various  neutralization  scenarios.  Under  Task  1.4,  a  particle  - 
turbulence  interaction  model  was  developed  for  application  to  the  spore  neutralization  model 
identified  under  Task  1.2.  Additionally,  in  coordination  with  Georgia  Tech,  the  spore 
neutralization  model  was  being  coupled  with  a  particulate  combustion  model  for  afterburning 
munitions.  Under  Task  1.5,  implementation  of  the  modeling  formulation  developed  under  Task 
1 .4  into  the  CRAFT  CFD®  flow  solver  was  initiated. 

1.3  Opportunities  for  Training  and  Professional  Development 

Nothing  to  report. 

1.4  Results  Dissemination  to  Communities  of  Interest 

Modeling  capabilities  of  the  CRAFT  CFD®  code  with  respect  to  biological  AD  turbulent  rate 
enhancement  developed  under  the  Basic  Research  program  were  highlighted  in  various  review 
meetings  of  other  DTRA  programs  in  which  CRAFT  Tech  was  involved. 
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2.0  YEAR  2  SUMMARY 


2.1  Project  Major  Goals 

The  primary  goals  for  the  Year  2  work  period  under  this  subeontraet  effort  were  as  follows: 

1 .  Develop  and  implement  a  turbulence  model  extension  to  account  for  turbulent  -  particle 
interactions  on  heterogeneous  combustion  of  afterburning  explosives  (ABX)  with  a 
particular  focus  on  aluminized  charges.  This  modeling  effort  focused  on  developing  a 
tractable,  production  level  capability  that  may  be  applied  within  the  context  of  routine, 
large  scale  analysis. 

2.  Integrate  modeling  formulations  for  spore  -  turbulence  interactions,  turbulence  -  particle 
interactions  for  heterogeneous  combustion,  and  hybrid  Reynolds  averaged  Navier- 
Stokes/large  eddy  simulation  (RANS/LES)  turbulence  modeling  for  momentum  and  heat 
transport  within  a  unified  computational  framework.  This  unified  framework  supports  a 
production  level  simulation  tool  specialized  for  blast  applications.  This  software  tool  will 
provide  analyst  with  the  capability  to  accurately  and  efficiently  assess  spore 
neutralization  scenarios  that  employ  ABX  charges.  This  capability  is  specialized  for 
computational  efficiency  within  the  context  of  routine,  large  scale  applications. 

3.  Develop  and  implement  turbulence  model  extensions  to  address  turbulence  -  chemistry 
interaction  for  gas  phase  combustion  within  the  context  of  production  level  applications. 
This  includes  a  two  phase  approach  to  upgrade  the  computational  framework  developed 
under  this  subcontract.  The  first  phase  focuses  on  implementation  of  a  first  order 
modeling  formulation  that  is  generally  applicable  to  account  for  the  primary  effect  of 
turbulent  combustion  on  gas  phase  chemical  production.  The  second  phase  of  this 
approach  focuses  on  investigating  the  application  of  a  more  advanced  modeling 
formulation  under  development  at  Georgia  Tech  within  production  level  applications. 
For  this  second  phase  effort  to  be  successful,  an  alternative  approach  must  be  found  to 
apply  this  modeling  formulation  in  a  more  efficient  and  computationally  tractable 
manner. 

2.2  Accomplishments 

During  the  Year  2  work  period,  efforts  focused  on  Tasks  1.5,  2.3  and  2.4.  Under  Task  1.5,  the 
development  and  implementation  of  a  particle  -  turbulence  interaction  model  for  spore 
neutralization  was  completed.  This  formulation  was  implemented  into  the  production  version  of 
the  CRAFT  CFD®  flow  solver  that  has  been  developed  for  application  to  blast  and  biological 
warfare  (BW)  neutralization  predictions.  Under  Task  2.3,  a  unified  modeling  formulation  was 
developed  and  implemented  to  fully  integrate  and  couple  the  modeling  formulations  developed 
under  the  previous  work  period’s  efforts.  This  effort  consisted  of  integrating  the  spore 
interaction  model  of  Task  1.4,  a  turbulence  -  particle  interaction  model  for  particulate 
combustion,  and  a  hybrid  Reynolds  averaged  Navier-Stokes/large  eddy  simulation  (RANS/LES) 
subgrid  model.  With  this  unified  formulation,  spore  neutralization  strategies  may  be  assessed  for 
afterburning  munitions  within  the  context  of  production  or  engineering  level  evaluations. 
Testing  of  this  unified  formulation  for  a  multi -room  blast/neutralization  problem  was  also  carried 
out.  A  key  finding  of  this  model  evaluation  effort  was  that  under  certain  conditions  turbulence  - 
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particle  interactions  significantly  reduce  spore  neutralization  rates.  This  model  prediction  has 
important  implications  for  the  assessment  of  prospective  neutralization  strategies.  Neglecting  the 
effect  of  turbulent  fluctuations  on  the  spore  destruction  rate  could  significantly  over  estimate  the 
effectiveness  of  a  neutralization  strategy.  As  a  result,  it  is  important  to  include  the  effect  of  these 
fluctuations  on  any  predictive  assessment  to  provide  a  more  conservative  estimate  of 
neutralization  effectiveness. 

Under  Task  2.4,  extensions  of  this  modeling  formulation  for  gas-phase  turbulence  chemistry 
interactions  were  initiated.  The  Eddy  Dissipation  Concept  (EDC)  formulation  for  turbulence  - 
chemistry  interactions  was  implemented  within  CRAFT  CFD®.  Also,  efforts  were  initiated  to 
investigate  a  more  advanced  combustion  modeling  formulation  based  on  the  linear-eddy  model 
(LEM)  in  coordination  with  Georgia  Tech.  Application  of  the  LEM  within  the  context  of  large 
scale  blast  applications  focused  on  development  of  a  tractable  formulation  based  of  stochastic 
model  parameterization.  Within  such  an  approach,  closure  statistic  generated  by  the  LEM  are 
parameterized  in  terms  of  a  reduced  set  of  variables  and  stored  within  a  model  database.  This 
database  may  then  be  deployed  within  a  flow  solver  to  account  for  turbulence  -  chemistry 
interaction  based  on  the  advanced  formulation  of  the  LEM,  but  at  a  small  fraction  of  the 
computational  cost. 

2.3  Opportunities  for  Training  and  Professional  Development 

Nothing  to  report. 

2.4  Results  Dissemination  to  Communities  of  Interest 

Modeling  capabilities  of  the  CRAFT  CFD®  code  with  respect  to  biological  AD  turbulent  rate 
enhancement  developed  under  the  Basic  Research  program  were  highlighted  in  various  review 
meetings  of  other  DTRA  programs  in  which  CRAFT  Tech  is  involved  as  well  as  at  a  C-WMD 
workshop  attended  by  the  broader  C-WMD  community  directly  supporting  DTRA.  Awareness 
on  the  potential  enhancing/inhibiting  effects  of  turbulent  mixing  on  a  number  of  finite-rate 
physical/chemical  processes  was  raised,  including  for  chemical  AD. 
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3.0  YEAR  3  SUMMARY 


3.1  Project  Major  Goals 

Under  this  subcontract  effort,  the  primary  goal  for  the  Year  3  work  period  was  to  perform 
verification  and  validation  (V&V)  of  the  predictive  capabilities  of  the  CRAFT  CFD®  code  with 
respect  to  realistic  blast  event/agent  defeat  (AD)  scenarios.  Specifically,  driven  by  DTRA's 
requirement  to  obtain  greater  confidence  in  the  accuracy  of  computational  simulations  involving 
advanced  high  energy  explosives  in  closed  compartments/structures  and  agent  defeat  simulations 
with  lower  energy  release  and  dispersion  for  neutralizers,  the  overarching  technical  objective 
was  to  assess  the  maturity  of  the  CRAFT  CFD®  code  in  capturing  the  effects  of  turbulent  mixing, 
disparate  time-scales,  multi-phase  physics  and  scenario  input  uncertainty  on  AD  in  a  complex 
geometry. 

3.2  Accomplishments 

During  the  Year  3  work  period,  efforts  focused  on  Tasks  3.3  and  3.4  pertaining  to  the 
verification  and  validation  (V&V)  of  the  modeling  capabilities  of  the  CRAFT  CFD®  code  with 
respect  to  agent  dispersion  and  AD  predictions,  respectively.  In  coordination  with  DTRA,  a 
specific  AD  test  configuration  of  interest  to  DTRA  was  selected.  Time  accurate  simulations  were 
conducted  of  a  realistic  AD  event  in  support  of  large  scale  test  activities.  As  part  of  this  analysis, 
comparison  against  measured  test  data  and  quantification  of  AD  effectiveness  of  an  enhanced 
high-energy  explosive  formulation  (compared  to  a  baseline  formulation)  were  carried  out. 
Details  of  this  analysis  were  provided  directly  to  DTRA  in  a  separate  technical  document. 

3.3  Opportunities  for  Training  and  Professional  Development 

Nothing  to  report. 

3.4  Results  Dissemination  to  Communities  of  Interest 

Modeling  capabilities  of  the  CRAFT  CFD®  code  with  respect  to  AD  turbulent  rate  enhancement 
developed  under  the  Basic  Research  program  have  been  highlighted  in  various  review  meetings 
of  other  DTRA  programs  in  which  CRAFT  Tech  is  involved.  Moreover,  awareness  on  the 
potential  enhancing/inhibiting  effects  of  turbulent  mixing  on  a  number  of  finite-rate 
physical/chemical  processes  was  raised  to  DTRA  and  the  broader  C-WMD  community.  With  the 
enhanced  fidelity  upgrades  and  V&V  of  the  CRAFT  CFD®  code  for  AD  scenarios  in  place, 
improved  predictive  capabilities  are  available  to  the  C-WMD  community  to  address  current 
threats  and  challenges,  and  design  enhanced  explosive  formulations.  Last,  opportunities  were 
identified  for  transitioning  high-fidelity  prediction  results  to  C-WMD  engineering  tools. 
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4.0  YEAR  4  SUMMARY 


4.1  Project  Major  Goals 

Driven  by  DTRA's  requirement  to  obtain  greater  confidence  in  the  accuracy  of  computational 
simulations  involving  advanced  high  energy  explosives  in  closed  compartments/structures  to 
support  agent  defeat  (AD)  with  lower  energy  release  and  dispersion  for  neutralizers,  the 
overarching  technical  objective  for  Year  4  of  this  program  was  to  assess  the  maturity  of  CRAFT 
Tech's  high-fidelity  CRAFT  CFD®  code  in  capturing  the  effects  of  turbulent  mixing,  disparate 
time-scales  and  multi-phase  physics  characteristic  of  the  AD  problem.  In  Year  4  of  this 
subcontract,  the  verification  and  validation  (V&V)  activities  of  the  CRAFT  CFD®  code  initiated 
in  Year  3  continued.  Particular  emphasis  was  placed  on  supporting  and  complementing  DTRA- 
funded  AD  testing  programs  and  also  in  addressing  AD  diagnostics  and  safety  requirements  for 
AD  tests  to  be  planned. 

4.2  Accomplishments 

During  the  Year  4  work  period,  efforts  continued  to  focus  on  Tasks  3.3  and  3.4  pertaining  to 
V&V  of  the  modeling  capabilities  of  the  CRAFT  CFD®  code  with  respect  to  AD  predictions.  In 
coordination  with  DTRA  and  leveraging  on  CRAFT  Tech's  active  involvement  in  the  AD 
Modeling  and  Simulation  (M&S)  community,  test  configurations  of  interest  to  DTRA  were 
evaluated  with  respect  to  V&V  as  well  as  with  respect  to  the  design  of  new  AD  tests. 
Specifically: 

•  In  preparation  of  chemical  AD  tests  with  simulants  in  the  small-scale  two-chamber  test 
configuration  at  a  Navy  facility,  a  parametric  study  of  the  relevant  design  variables  (e.g., 
amounts  of  HE,  additives,  simulant,  etc.)  was  carried  out  to  narrow  down  the  parameter 
ranges  to  ensure  safety  and  attain  a  representative  AD  scenario  at  the  laboratory- scale. 
Based  on  these  findings  and  on  discussions  with  Navy  personnel,  a  preliminary  test 
design  was  identified,  for  which  a  high-fidelity  simulation  with  the  CRAFT  CFD®  code 
was  carried  out.  Key  aspects  of  interest  in  this  simulations  were  (i)  the  characterization  of 
chemical  agent  simulant  neutralization  and  quantification  of  amounts  released  from  the 
blast  chamber  into  the  exhaust  chamber,  and  (ii)  characterization  of  pressure  and 
temperature  levels  in  both  chambers.  Numerical  results  have  been  delivered  to  the  Navy 
for  comparison  with  their  experimental  data.  Details  of  this  high-fidelity  CFD  analysis 
were  documented  and  were  provided  directly  to  DTRA  in  a  separate  technical  document. 

•  As  a  member  of  the  AD  M&S  community,  CRAFT  Tech  is  involved  in  high-fidelity 
V&V  activities  and  AD  test  planning  support.  For  a  geometry  of  interest,  additional  side 
calculations  were  performed  to  exercise  features  or  capabilities  of  the  CRAFT  CFD® 
code  pertinent  to  turbulent  mixing  and  AD  effectiveness: 

o  Since  the  temperature  sensed  by  the  agent  as  it  is  dispersed  by  turbulent  mixing  is 
a  critical  AD  parameter,  a  diagnostic  approach  under  development  and  stemming 
from  DTRA  Basic  Research  program  is  represented  by  micron-sized 
thermometers.  As  a  demonstration  of  the  ability  of  high-fidelity  modeling  to 
provide  AD  test  planning  support  in  estimating  the  temperature  history 
experienced  by  these  micro-thermometers  and  helping  with  the  interpretation  of 
the  micro-thermometer  measurement  data,  a  computational  particle  tracking 
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exercise  was  performed  in  a  post-detonation  environment  subject  to  turbulent 
mixing  to  provide  time-resolved  temperature  histories  and  particle  trajectories  of 
selected  tracers. 

o  A  second  computational  exercise  was  performed  with  regard  to  deploying  a 
simple,  computationally-tractable  turbulent  combustion  model  like  the  Eddy 
Dissipation  Concept  (EDC  -  refer  to  Year  2  activities)  in  the  same  geometry  to 
highlight  agent  neutralization  rate  sensitivities. 

4.3  Opportunities  for  Training  and  Professional  Development 

Nothing  to  report. 

4.4  Results  Dissemination  to  Communities  of  Interest 

Modeling  capabilities  of  the  CRAFT  CFD®  code  with  respect  to  AD  turbulent  rate  enhancement 
developed  under  the  Basic  Research  program  have  been  highlighted  in  various  review  meetings 
of  other  DTRA  programs  in  which  CRAFT  Tech  is  involved.  Moreover,  awareness  on  the 
potential  enhancing/inhibiting  effects  of  turbulent  mixing  on  a  number  of  finite-rate 
physical/chemical  processes  was  raised  to  DTRA  and  the  broader  C-WMD  community.  With  the 
enhanced  fidelity  upgrades  and  V&V  of  the  CRAFT  CFD®  code  for  AD  scenarios  in  place, 
improved  predictive  capabilities  are  available  to  the  C-WMD  community  to  address  current 
threats  and  challenges,  and  design  enhanced  explosive  formulations.  Also,  opportunities  were 
identified  for  transitioning  high-fidelity  prediction  results  to  C-WMD  engineering  tools.  Fast, 
CRAFT  Tech  has  actively  participated  in  the  December  2015  Simulant  Chemistry  Workshop 
organized  by  DTRA.  This  was  an  excellent  opportunity  to  further  strengthen  the  link  between 
DTRA's  Basic  Research  portfolio  and  the  applications  side  of  DTRA's  activities. 
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5.0  YEAR  1  ACTIVITIES 


5.1  Task  1.2:  Database  Generation  of  BW  Agent  Neutralization  Models 

5.1.1  Overview 

In  this  task,  a  survey  of  the  state-of-the-art  in  biological  warfare  (BW)  agent  defeat  (AD) 
modeling  was  performed.  Relevant  BW  agent  defeat  models  for  thermal  and/or  chemical 
neutralization  were  selected  and  included  in  a  spore  kill  database  for  use  within  CRAFT  Tech’s 
high-fidelity  multi-phase  interaction  framework.  It  was  found  that  for  large-scale  simulations  of 
blast  events,  global  models  that  rely  on  empirical  correlations  are  suitable.  The  spore  kill 
database  targets  primarily  neutralization  of  Bacillus  anthracis  {B.a.)  and  Bacillus  Thuringiensis 
(B.t.)  exposed  to  a  high  temperature  gas  environment  (HTGE).  Experimental  measurements  on 
their  rate  of  demise  are  generally  performed  under  controlled  conditions,  uniform  in  space  and 
constant  in  time,  and  data  is  then  correlated  in  the  form  of  an  Arrhenius-like  neutralization  rate 
by  assuming  zero-order  or  first-order  kinetics,  as  done  in  chemical  kinetics  modeling.  During 
agent  defeat  events,  e.g.,  a  HE  blast  within  a  hardened  target,  conditions  are  highly  variable. 
Therefore,  within  a  modeling  framework,  the  form  of  the  neutralization  rate  should  be 
sufficiently  generalized  to  account  for  dependency  on  various  parameters,  e.g.,  humidity  level, 
through  a  global  representation.  Rate  data  can  be  adjusted  for  specific  conditions  or  replaced  as 
new  data  or  improved  models  become  available. 

5.1.2  Formulation 

HTGE  or  a  corrosive  agent  impact  spore  viability,  i.e.,  their  ability  to  form  colonies.  Spore 
neutralization  is  generally  quantified  in  terms  of  log  10  reduction  in  viability.  For  instance,  if  no 
is  the  initial  spore  number  density  number  and  n  is  the  spore  number  density  after  exposure  to 
HGTE,  the  viability  ratio  is  defined  as: 


V  =10“''  (1) 

/% 

While  a  five-decade  viability  reduction  is  generally  accepted  (r  =  5),  a  kill  effectiveness  of  r  =  9 
is  desirable.  Assuming  a  first-order  global  representation  of  neutralization,  the  evolution  in  time 
of  the  viability  ratio  is  governed  by  the  ordinary  differential  equation: 


where  the  neutralization  rate  k  is  expressed  in  an  Arrhenius-like  form: 


(2) 


k  =  aT'"e  ^  (3) 

The  rate  parameters  a,  corresponding  to  a  collision  frequency,  P,  corresponding  to  an  activation 
energy,  and,  m,  representing  a  temperature  exponent  are  determined  from  experiments. 
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5.1.3  Literature  Survey 

Past  measurements  of  spore  thermal  inactivation  have  focused  on  relatively  long  exposure  times, 
spanning  from  minutes  to  several  hours,  in  moderate  HGTE.  Starting  from  the  early  work  of 
Femelius  on  B.a.  and  Bacillus  Subtilis  (B.s.)  [1]  to  the  more  recent  comprehensive  survey  of 
Spotts  Whitney  at  CDC  on  B.a.  [2],  the  temperature  range  considered  remains  low.  Table  1 
summarizes  the  data  collected  for  moist  heat  (90  C  -  120  C)  and  for  dry  heat  (140  C  -  200  C). 


Table  1.  Summary  of  data  surveyed  by  Spott  Whitney  [2]. 


Temperature 

Time 

Inoculum  size 

Temperature 

Time 

Inoculum  size 

Boiling 

Dry  heat 

100"  c 

10  min 

3x  10® 

140"  C 

>90  min 

6x  10®  to  1.2x  10‘ 

5  min 

7.5x  10® 

150"  C 

10  min 

6x  10®  to  1.2x  10‘ 

Moist  heat 

160"  C 

10  min 

6x  10®  to  1.2x  10‘ 

90"  C 

20  min 

1.2x  10® 

180"  C 

2  min 

6x  10®  to  1.2x  10‘ 

90"  Cto91"  C 

60  min 

3x  10® 

190"  C 

1  min 

6x  10®  to  1.2x  10‘ 

100"  C 

10  min 

1.2x  10® 

200"  C 

30  sec 

6x  10®  to  1.2x  10‘ 

100"  CtolOI"  C 

17  min 

lx  10® 

105"  C 

10  min 

3x  10® 

120"  C 

15  min 

2.4x  10® 

For  HGTE  applications  generated  by  a  HE  blast,  e.g.,  munitions  explosion,  characteristic  spore 
exposure  time  scales  are  less  than  1  s  and  the  peak  temperatures  reached  span  from  several 
hundreds  to  more  than  a  thousand  degrees  K.  Figure  1  illustrates  the  distribution  of  HGTE 
temperatures  and  heating  times  encountered  in  HE  applications  and  in  long  duration  events  (from 
Burggraf  at  the  Air  Force  Institute  of  Technology  (AFIT)  [3] [4]).  Research  programs  at  AFIT 
have  focused  on  short  duration  thermal  inactivation,  both  experimental  characterization  and 
numerical  modeling  accounting  for  the  spore  structure  [3]  [4],  Particular  emphasis  was  given  to 
hydrolysis  reactions  and  effects  on  the  spore  DNA. 


o 


B.a.  Spotts-Whitney,  et  al. 
B.a.  Fernellus  et  al. 

B.a.  Goetz 
B.a.  hlawkins 
B.t.  hlawkins 
B.t.  Battell  5  decades 
B.t.  Battell  2  decades 
HE  Region 
Incendiary  Region 


Figure  1:  Summary  of  distribution  of  temperature  and  thermal  exposure  duration  (from 

Burggraf,  AFIT). 
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Similarly,  Levin  et  al.  conducted  detailed  modeling  of  B.a.  spores  in  HGTE  to  support  the 
interpretation  of  experimental  findings  [5].  Specifically,  CFD  and  DSMC  simulations  were 
carried  out  to  analyze  the  transient  heat  transfer  to  an  individual  spore  based  on  the  estimation  of 
the  heat  transfer  coefficient.  Disagreement  between  gas-surface  heat  conduction  results  and 
experiments  (the  latter  indicated  a  much  longer  deactivation  time  scale)  prompted  to  speculate 
the  effect  of  internal  pressure  exerted  on  the  spore  walls  due  to  heating  or  the  clumping  of  the 
spores,  thereby  providing  a  shielding  effect  and  reducing  the  heat  transfer  rate  to  the  spore. 
Although  such  detailed  modeling  analysis  accounting  for  the  spore  structure  and  clump 
configuration  are  generally  not  usable  directly  in  large-scale  simulations  of  blast  events,  they  are 
very  valuable  in  the  development  of  reduced-order  models.  In  the  future,  these  reduced-order 
models,  consistent  with  a  generalized  neutralization  rate  formulation,  can  be  integrated  into 
CRAFT  Tech's  modeling  framework. 

Empirical  correlations  for  B.t.  demise  have  been  formulated  in  the  423  K  -  1210  K  temperature 
range  at  Battelle  [6]  [7].  Assuming  a  five  decade  viability  reduction  for  neutralization,  the 
residence  time  x  in  HTGE  required  as  a  function  of  the  air  temperature  T  was  derived.  Similarly, 
correlations  for  the  concurrent  effect  of  HTGE  and  corrosive  gases  have  been  defined  in  terms  of 
a  single  Arrhenius  kinetics  relations  [8].  The  primary  argument  is  that  HTGE  directly  attacks  the 
spores  and  at  the  same  time  the  high  temperature  accelerates  the  chemical  attack  of  the  corrosive 
agents.  Under  these  conditions  with  biocidal  material,  reaction  orders  may  not  be  purely  zero-th 
or  first  order.  A  combined  neutralization  rate  formula  for  Thiokol  TP-H1246  propellant  was 
derived. 


More  recently,  SwRI  has  performed  extensive  experimental  work  in  support  of  further 
development  of  the  Empirical  Lethality  Methodology  (ELM)  [35].  The  Air  Force  Nuclear 
Weapons  Center,  Capabilities  and  Integration  Directorate  (AFNWC/XR)  maintains  an  agent 
defeat  database  addressing  thermal,  chemical  and  radiative  neutralization.  Testing  was  performed 
on  Bt  var.  kurstaki,  with  particular  emphasis  on  thermal  neutralization  and  on  effect  of  spores  in 
hydrated  form,  for  instance  moist  spores  or  spores  in  a  suspension.  A  test  matrix  was  constructed 
targeting  environment  temperatures  of  200  C,  300  C  and  400  C,  and  relative  humidity  (RH) 
levels  of  27%,  48%,  82%  and  96%.  A  substantial  amount  of  data  was  collected  based  on 
repeated  testing  attempting  to  correlate  key  variables,  namely  (i)  environment  temperature  T,  (ii) 
RH  and  (iii)  residence  time  x  of  the  spores  in  the  HTGE  to  the  viability  ratio  r,  characterized  in 
terms  of  Blood  Agar  Plate  (BAP)  log  10  kill.  The  raw  data  was  processed  into  an  apparent 
Arrhenius-like  rate  assuming  a  first-order  reaction  rate  order  for  the  rate  of  spore  demise: 


Ink^ 


-6  —  +  \na 
T 


(4) 


Specifically,  a  linear  trend  is  expected  between  the  inverse  of  the  environment  temperature  and 
the  logarithm  of  the  neutralization  rate.  Here,  the  neutralization  rate  kRH  is  computed  from  the 
exposure  time  x  and  viability  ratio  r  based  on  the  analytical  solution  to  the  neutralization 
equation: 


kR„iT)  = 


rlnlO 


(5) 


T 

Despite  the  significant  scatter  observed,  several  sets  of  Arrhenius-like  rate  parameters  were 
identified  at  various  levels  of  RH.  This  rate  data  for  kRH  in  terms  of  a  and  p  can  be  deployed 
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directly  into  CRAFT  Tech's  high-fidelity  modeling  framework  and  has  been  incorporated  into 
the  spore  kill  database. 

5.1.4  Additional  Sensitive  Sources  for  Spore  Neutralization  Rates 

Additional  data  is  potentially  available  from  Government-generated  agent  defeat  databases.  For 
instance,  these  may  include  the  Empirical  Lethality  Methodology  (ELM)  database  [35],  the  US 
Army’s  Edgewood  Chemical  Biological  Center  (ECBC)  Agent/Simulant  Knowledge  (ASK) 
database  [10],  and  the  data  generated  from  an  SBIR  program  with  Air  Force  in  2003-2004  [11]. 
Although  this  data  was  not  evaluated  as  part  of  the  present  task,  it  is  expected  to  be  compatible 
with  the  generalized  spore  neutralization  form  adopted  in  CRAFT  Tech's  modeling  framework. 

5.1.5  Additional  Physico-Chemical  Properties 

Moreover,  additional  physico-chemical  properties  beyond  neutralization  rates  are  required  for  an 
accurate  characterization  of  spore  behavior  within  a  high-fidelity  modeling  framework.  For 
instance,  these  properties  entail  spore  density,  specific  heat,  size  distribution,  specific  volume. 
For  B.a.  spores.  Levin  et  al.  have  estimated  density  at  about  1300.0  kgW  and  specific  heat 
capacity  at  2500.0  J/kg/K  [5].  The  U.S.  Army's  ECBC  has  characterized  the  size  distribution  of 
various  strains  of  B.a.  spores  as  well  as  several  other  spores,  including,  B.t,  Bacillus  subtilis. 
Bacillus  cereus,  etc.  [12].  For  each  strain  and  spore,  they  provide  mean  length,  diameter,  aspect 
ratio  and  volume  as  well  as  the  range  of  variability.  On  average,  a  B.a.  spore  has  a  length  of 
slightly  above  1  pm  and  a  diameter  of  above  0.8  pm. 

5.2  Task  1.4:  Turbulence  Model  Extensions  for  Biological  Agent  Neutralization 

5.2.1  Particle  -  Turbulence  Interaction  Modeling  for  BW  Agent  Neutralization 

In  cooperation  with  Georgia  Tech,  CRAFT  Tech  has  been  focused  on  developing  a  particle  - 
turbulence  interaction  model  for  biological  agent  neutralization  for  application  to  large  scale 
problems  of  interest  to  DTRA.  One  primary  requirement  for  this  formulation  is  that  it  be 
tractable  for  routine  application  to  realist,  large  scale  problems.  With  this  requirement  in  mind,  a 
parameterized  approach  was  chosen  for  the  development  of  the  formulation.  Within  this 
approach,  the  statistics  from  the  physical  models,  such  as  those  under  development  at  Georgia 
Tech,  are  parameterized  into  an  efficient  database.  In  this  manner,  a  large  scale  computation  will 
access  the  model  database  to  obtain  the  required  statistics  in  an  efficient  and  cost  effective 
manner.  This  approach  is  in  contrast  to  directly  employing  the  physical  model  within  the 
computation,  which  will  be  intractable  for  large  scale  problems  of  DTRA  interest. 

Within  this  overall  approach,  the  modeling  efforts  thus  far  have  focused  on  two  aspects  of 
particle  -  turbulence  interactions.  First,  a  parameterized  model  for  spore  neutralization  has  been 
developed  to  account  for  turbulence  flow  fluctuations  on  the  neutralization  rate.  Second, 
following  parallel  efforts  at  Georgia  Tech,  a  separate  particle  -  turbulence  interaction  model  has 
been  developed  for  aluminum  particle  combustion  for  application  to  afterburning  munitions  that 
could  be  used  for  spore  neutralization.  These  two  modeling  formulations  will  be  coupled 
together  to  extend  the  applicability  of  the  neutralization  modeling  for  application  to  afterburning 
munitions. 
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The  following  subsections  of  this  report  summarize  the  modeling  formulations  for  both  the  spore 
neutralization  and  the  aluminum  combustion. 

5.2.2  Spore  Particle  Interaction  Modeling 

As  discussed  earlier,  correlations  for  the  rate  of  spore  neutralization  may  be  specified  in  the 
general  form, 


w,=-pY,aT’”exp(-l5/T)  (6) 

where  p  is  the  gas  phase  density,  7^  is  the  live  spore  mass  concentration,  T  is  the  gas  phase 
temperature,  and  m,  a  and  (3  are  coefficients  specified  from  experimental  correlations.  In  this 
equation,  the  spore  temperature  is  assumed  to  be  in  equilibrium  with  the  gas  phase  so  that  the  gas 
phase  temperature  appears  in  Eqn.  (24).  Application  of  Eqn.  (24)  to  turbulent  flow  requires 
characterizing  the  effects  of  spore  concentration  and  gas  phase  property  fluctuations  on  the  mean 
neutralization  rate  which  may  be  written  as, 

^  /  T )P(p,  ,  T )dpdY^dT  (7) 

where  P  is  the  joint  probability  density  function  of  the  gas  phase  density  and  temperature,  and 
the  live  spore  mass  concentration.  Given  P,  the  Eqn.  (25)  may  be  evaluated  to  characterize  the 
effect  of  turbulent  fluctuations  on  the  mean  neutralization  rate. 

The  pdf  in  Eqn.  (25)  may  be  specified  by  either  direct  evaluation  or  a  parameterization  strategy. 
The  direct  evaluation  approach  seeks  to  predict  the  shape  of  the  pdf  through  the  solution  of  a 
transport  equation  for  P  based  on  fundamental  principles.  This  approach  is  comprehensive,  but 
intractable  for  all  but  a  limited  range  of  simplified  problems.  Alternatively,  the  effect  of 
turbulent  fluctuations  within  Eqn.  (25)  may  be  characterized  through  a  parameterization  of  P  in 
terms  of  a  reduced  set  of  variables.  Through  such  an  approach,  a  tractable  modeling  formulation 
may  be  developed  that  may  be  applied  to  large  scale  problems  of  interest.  This  type  of 
parameterization  approach  was  chosen  to  characterize  the  effect  of  turbulent  fluctuations  on  the 
spore  neutralization  rate. 

This  work  period,  a  parameterized  modeling  approach  was  developed  to  model  Eqn.  (25). 
Within  this  approach,  the  pdf  P  was  parameterized  using  an  assumed  pdf  formulation  [19] 
assuming  statistical  independence  between  the  density,  temperature  and  spore  concentration. 
With  this  assumption,  the  pdf  in  Eqn.  (25)  may  be  written  as. 


P(p.T.Y,)~PJp-p)Pr(T-T)P,(Y,  -7J  (8) 

where  P^,  Pj.,  and  Py  are  the  marginal  pdfs  the  density,  temperature  and  spore  concentration. 

These  marginal  pdfs  are  then  specified  as  follows.  For  density,  a  delta  function  is  specified  for 
Pp  [19].  For  Pj.  and  Py ,  a  beta  function  [19]  is  used  which  is  parameterized  in  terms  of  it  first 

two  moments.  For  example,  using  a  beta  function,  Py  may  be  expressed  as. 


PyiY,)  = 


r(A)r(A) 


r(A  +7^2) 


(9) 
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with, 


-1 


,  A  =(1-^1) 


Y,(\-Y,) 


-1 


-  '2 


(10) 


where  T (x)  is  the  gamma  function.  As  seen  in  these  equations,  the  beta  function  for  has  been 


parameterized  in  terms  of  the  mean  spore  concentration  (  7^  )  and  its  variance  A  similar 

expression  may  be  developed  for  Pj.  that  is  parameterized  in  terms  of  the  mean  temperature  and 
temperature  variance. 


With  this  parameterization  of  Eqn.  (26),  the  mean  spore  neutralization  rate  in  Eqn.  (25)  may  be 
evaluated.  With  this  formulation,  a  pdf  table  generator  code  was  developed  this  work  period  to 
generate  a  database  for  the  mean  spore  rate  as  a  function  of  the  mean  temperature,  temperature 
variance,  mean  spore  concentration,  and  spore  concentration  variance.  This  database  may  then 
be  used  to  specify  the  mean  spore  rate  (Eqn.  (25))  in  the  spore  transport  equation. 


This  parameterized  formulation  additionally  requires  the  specification  of  the  temperature  and 
spore  concentration  variances  that  appears  in  the  beta  representations  of  Pj.  and  Py .  These 

quantities  are  specified  from  addition  transport  equations.  For  the  temperature  variance,  a 
variable  turbulent  Prandtl  number  model  [20]  that  is  applicable  to  reacting  flows  is  used.  For  the 
spore  fluctuations,  a  spore  concentration  variance  equation  was  developed  that  is  similar  to  the 
species  concentration  variance  methodology  of  Gaffney,  et  al.  [21].  This  equation  includes  an 

additional  term  of  the  form  7^w^ ,  which  is  also  modeling  using  the  assumed  pdf  methodology  of 

Eqn.  (26).  The  parameterized  model  of  this  term  is  also  stored  within  the  model  database  that  is 
generated  by  the  pdf  table  generator  that  was  developed  this  work  period. 

The  development  of  the  pdf  table  generator  for  the  mean  spore  rate  was  completed  this  work 
period.  Also,  this  work  period,  the  implementation  of  the  mean  spore  concentration  and 
concentration  variance  equations  into  the  CRAFT  CFD®  code  was  initiated  under  Task  1.5.  This 
formulation  is  being  implemented  into  the  version  of  the  CRAFT  CFD®  that  includes  the 
aluminum  combustion  model  for  afterburning  munitions.  This  modeling  formulation  is 
described  in  the  next  subsection. 


5.2.3  Particle  -  Turbulence  Interaction  Modeling  for  Aluminum  Particulates 

For  turbulent  flow,  velocity  and  composition  fluctuations  will  influence  the  burning  rate  of  solid 
aluminum  particulates  within  afterburning  munitions  that  may  be  used  for  spore  neutralization. 
To  account  for  these  fluctuations,  two  approaches  may  be  considered.  The  first  is  a  microscopic 
approach  that  seeks  to  model  the  direct  interaction  of  the  flow  with  the  flame  structure 
surrounding  the  particles  from  a  first  principles  perspective.  The  second  is  a  macroscopic 
approach  that  seeks  to  account  for  flowfield  fluctuations  on  the  modeled  laminar  burning  rate. 
The  microscopic  approach  is  more  fundamental,  but  intractable  from  the  perspective  of  applying 
it  routinely  to  problems  of  practical  interest.  The  macroscopic  approach  is  more  pragmatic  and 
empirical,  but  provides  a  viable  alternative  to  develop  a  tractable  formulation.  As  a 
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consequence,  modeling  efforts  focused  on  the  development  of  a  macroseopie  particle  - 
turbulenee  interaction  model  that  may  be  routinely  applied  to  problems  of  interest. 

For  the  macroseopie  approach,  the  model  sought  to  aecount  for  the  effects  of  turbulent 
fluetuations  on  the  modeled  laminar  particle  burning  rate.  Within  CRAFT  Teeh’s  Eulerian 
dispersed  phase  particle  eombustion  modeling  formulation,  the  laminar  aluminum  particle  mass 
density  eonsumption  rate  is  expressed  as. 


mp=-2jrrlNp^^  (H) 

at 

where  Pp  is  the  particle  mass  density,  pm  is  the  material  density,  N  is  the  particle  number  density, 
and  Vp  and  dp  are  the  partiele  radius  and  diameter,  respeetively.  The  time  derivative  in  Eqn.  (30) 
is  also  expressed  as. 


_ 1 _ 

dt  ndpA^+md'^-^A^ 

with. 


(12) 


A,exp(-E,/RT^)X^j^\-8’^ 


(13) 


A^  = 


pOA  ^ 


eff 


1-f" 


(14) 


where  n  =  0.3,  m  =  1.8,  Ai  =  S.SxlO"*,  A2  =  7.35x10'*’,  e  =  0.05,  Eb  =  73.6  kJ/mol,  R  is  the  gas 
constant,  Tp  is  the  partiele  surfaee  temperature,  P  is  the  pressure,  T  is  the  gas  temperature.  The 
variable  Xef/is  the  sum  of  gas  phase  mole  fraetions  of  O2,  H2O  and  CO2  as  given  by. 


X^^f  =  +  0.6X^^^  +  0.22X^^^  (15) 

Eqns.  (30)  -  (34)  were  developed  from  the  Yetter’s  hybrid  eombustion  model  [22]  for  aluminum 
partiele  burning.  This  model  formulation  aceounts  for  both  the  diffusion-controlled  and 
kinetieally-eontrolled  burning  limits,  and  a  transition  between  the  two.  With  this  formulation, 
the  partiele  mass  density  eonsumption  rate  in  Eqn.  (30)  may  be  written  as  a  function  of  the 
partiele  and  gas  phase  properties  as. 


m,-f(r,.T,.P.T.X,)  (16) 

where  Xt  are  the  speeies  mole  fractions.  This  equation  is  a  representation  of  the  laminar  burning 
rate  model  within  CRAFT  Teeh’s  dispersed  phase  combustion  model  for  aluminum  particles. 

With  the  laminar  rate  model  given  by  Eqn.  (35),  a  macroseopie  level  model  for  the  effeet  of 
turbulent  fluetuations  on  the  rate  expression  could  be  represented  by, 

^ -Jpdf(r^.T^,P,T,Xt)f(r^.T,.P.T.Xt)dr/T/PdTdX^  (17) 

where  pdf  ( rp,Tp,P,T ,Xi^)  is  the  joint  probability  density  function  of  the  particle  and  gas  phase 
properties  and  the  overbar  represents  Reynolds  averaging.  With  Eqn.  (36),  a  parameterized 
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model  for  this  mean  rate  expression  has  been  developed.  As  discussed  earlier,  a  parameterized 
modeling  formulation  seeks  to  represent  the  model  statistics  in  terms  of  a  reduced  set  of 
variables.  These  model  statistics  may  then  be  stored  within  a  database.  A  flow  solver  may  then 
retrieve  the  model  statistics  from  the  database  in  a  fast  and  efficient  manner. 

To  develop  a  parameterized  model  for  Eqn.  (36),  it  should  first  be  recognized  that  the  particle 
and  gas  phase  properties  for  dispersed  phase  combustion  are  effectively  uncorrelated.  That  is, 
specific  gas  phase  properties  do  not  necessary  correspond  to  specific  particle  phase  properties. 
For  example,  high  instantaneous  values  of  gas  phase  temperature,  T,  do  not  necessarily 
correspond  to  high  instantaneous  values  of  particle  temperature,  or  specific  values  of  particle 
radius.  The  particle  and  gas  phase  properties  will  only  be  strongly  correlated  when  the  two 
phases  are  in  equilibrium.  As  a  result,  the  particle  and  gas  phase  properties  may  be  assumed  to  be 
statistically  independent.  With  this  assumption,  the  joint  pdf  in  Eqn.  (36)  may  be  expressed  as  a 
combination  of  the  marginal  pdfs  of  these  variables  as, 

pdf(r^.T^.P.T.X,)-pdf(r^.T^)pdf(P.T.X,)  (18) 

With  this  representation  of  the  joint  pdf,  an  assumed  pdf  model  [19]  may  be  used  to  develop  an 
approximation  for  the  RHS  of  Eqn.  (37).  To  develop  a  first  order  model,  the  pdf  of  the  particle 
properties  is  assumed  to  be  composed  of  delta  functions  at  the  mean  particle  values.  For  the  gas 
phase  pdf,  the  pressure  is  also  assumed  to  be  statistically  independent  of  the  temperature  and 
species.  These  assumptions  are  reasonable  for  the  present  application  because  the  particle  and 
gas  phase  properties  are  effectively  uncorrelated  as  discuss  earlier.  With  these  assumptions,  Eqn. 
(37)  may  be  rewritten  as, 


pdf(r^.T^.P.T.X,).d(r^ -r^)d(T^ -T^)pdf{P)pdf(T.X,)  (19) 

For  the  gas  phase  pdfs,  assumed  pdf  methods  typically  employ  a  delta  function  for  the  marginal 
pdf  of  pressure.  The  task  then  becomes  to  develop  a  parameterized  model  for  pdf  (T,Xi^) .  For 

afterburning  munitions,  particle  -  turbulence  interactions  will  primarily  occur  in  the  shear  layer 
region  between  the  ambient  fluid  and  the  blast  constituents  where  mixing  and  combustion  may 
occur.  Under  these  conditions,  the  mixture  gas  phase  properties  could  be  parameterized  in  terms 
of  mixture  fraction  (Z)  and  scalar  dissipation  (x)->  similar  to  the  parameterization  used  by 
Sankaran,  et  al.  [23]  for  gas  phase  combustion  and  for  flamelet  models  [19].  Consequently,  the 
local  temperature  and  species  within  the  shear  layer  may  be  written  as  T  =  T(Z,x,P)  and 

=  X,^(Z,x,P),  with  the  pressure  being  included  due  to  compressibility  effects  within  the 

flow.  This  parameterization  of  the  gas  phase  properties  may  be  accomplished  using  a  flamelet 
model  [19]  given  representative  thermodynamic  properties  of  the  flow  (i.e.,  fuel  and  oxidizer 
composition,  and  pressure  range).  With  this  representation,  and  the  aforementioned  assumption 
regarding  the  statistical  independence  of  the  pressure  and  temperature  and  species,  Eqn.  (38) 
then  becomes, 

pdf(r^  .T^.P.T.X,)~d(r^-7^)d(T^-%)d(P-P)S(x- x)pdf(Z)  (20) 

where  the  typical  assumptions  regarding  scalar  dissipation  from  assumed  pdf  methods  have  been 
employed  [19]. 
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Within  Eqn.  (39),  the  pdf  of  mixture  fraction  may  be  approximated  using  a  beta  function  [19]. 
The  mean  particle  density  consumption  rate  in  Eqn.  (39)  may  be  evaluated  as  a  function  of  the 
mean  flow  and  particle  properties  as, 

^,-g(r,.T^)h(P.X.Z.^)  (21) 

where  the  functions  g  and  h  are  independent  due  to  the  previous  assumption  that  the  particle  and 
flow  properties  are  uncorrelated.  Alternatively,  a  function  F  may  be  defined  as  the  ratio  of 
turbulent-to-laminar  burning  rate  as. 


F(r,T,P,x,Z,Z'\T)  = 


m 


(22) 


p,lam 


where  is  the  consumption  rate  defined  by  Eqn.  (35)  evaluated  based  on  mean  flow 

quantities  only.  The  term  is  the  value  of  the  consumption  rate  that  neglects  all  turbulent 

fluctuations,  or  that  is  the  effective  laminar  rate. 


With  Eqns.  (39)-(41),  a  database  for  the  function  F  may  be  evaluated  given  the  representative 
thermodynamics  conditions  of  the  problem.  This  database  may  then  be  deployed  within  a  CFD 
flow  solver  so  that  the  mean  particle  density  consumption  rate  may  be  evaluated  as, 


’n,-m^,„F(r^.T^.P.X.Z.Z'\T)  (23) 

By  defining  the  mean  consumption  rate  in  this  manner,  the  first  order  effect  of  turbulent 
flowfield  fluctuations  may  be  included.  To  apply  this  model,  the  CFD  flow  solver  must 

additionally  solve  standard  transport  equations  for  Z  and  Z'^ .  These  transport  equation  may  be 
specified  using  the  compressible  flow,  variable  turbulent  Schmidt  number  formulation  of 
Brinckman,  et  al.  [24]. 

With  this  modeling  formulation,  the  CFD  flow  solver  is  only  required  to  include  just  two 
additional  transport  equations.  The  multiplying  function  F  is  also  retrieved  from  a  database  in  a 
fast  and  efficient  manner.  As  a  result,  this  first  order  turbulence  -  particle  interaction  model  is 
computationally  inexpensive  to  employ  for  large  scale  applications.  The  effects  of  this  modeling 
formulation  on  aluminum  particle  combustion  may  be  seen,  for  example,  in  Figure  2  for  an 
unsteady,  high  speed  shear  layer.  This  figure  presents  the  reaction  enhancement  factor  F  in  Eqn. 
(41).  This  figure  shows  contours  of  this  function  in  the  ignition  region  of  the  aluminum 
particulates. 
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Figure  2:  Instantaneous  particle  reaction  enhancement  factor  contours  for  the  unsteady 
particle  burning  a  P  =  1  atm,  Tfud  =  2000  K  and  1  micron  particles 


This  aluminum  combustion  modeling  formulation  is  currently  being  coupled  with  the  spore 
combustion  model  described  in  the  previous  subsection  under  Task  1.5. 

5.3  Task  1.5:  Combined  Effort  for  Model  Integration  in  CRAFT  Tech  Codes 

Under  this  task,  efforts  were  initiated  to  integrate  the  aforementioned  spore  neutralization  model 
into  CRAFT  CFD®.  The  spore  modeling  formulation  is  being  integrated  within  the  aluminum 
combustion  model  that  has  been  developed  for  afterburning  munitions.  With  the  completion  of 
the  modeling  formulation  and  the  model  database  generator  under  Task  1.4,  implementation  of 
the  spore  transport  and  spore  concentration  variance  equations  into  CRAFT  CFD®  was  initiated. 
With  the  completion  of  this  implementation,  testing,  evaluation  and  validation  of  the  formulation 
will  be  carried  in  Year  2. 
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6.0  YEAR  2  ACTIVITIES 


6.1  Task  1.5:  Combined  Effort  for  Model  Integration  in  CRAFT  Tech  Codes 

Under  this  task,  the  development  and  implementation  of  the  particle  -  turbulence  interaction 
model  for  spore  neutralization  that  was  developed  in  Year  1  under  Task  1.4  was  completed. 
This  formulation  characterizes  the  influence  of  turbulent  flow  scalar  fluctuations  on  the  evolution 
of  spore  particulates  and  the  spore  neutralization  rate.  This  model  employs  the  neutralization 
rate  formulation  identified  under  Task  1.2.  The  effect  of  turbulent  flow  fluctuations  on  the 
evolution  of  the  spores  is  included  within  a  computationally  tractable  parameterization  approach. 
Within  such  an  approach,  the  effects  of  scalar  fluctuations  on  the  mean  or  subgrid  averaged  spore 
neutralization  rate  is  included  by  parameterizing  spore  rate  fluctuation  statistics  from  a 
computational  sub-model.  These  statistics  are  then  organized  within  a  database  which  is  used  for 
a  run-time  calculation.  In  such  a  manner,  the  model  statistics  required  by  the  CFD  flow  solver 
are  obtained  from  the  model  database  and  are  not  computed  “on  the  fly”.  As  such,  this 
formulation  is  computationally  tractable  for  large  scale,  production  level  applications  because  the 
computational  cost  of  accessing  the  database  is  much  less  than  the  direct  application  of  the  sub¬ 
model  at  run-time.  This  spore  particle  -  turbulence  interaction  modeling  formulation  is 
summarized  in  a  following  sub-section  of  this  report. 

During  the  Year  2  work  period,  the  development  and  implementation  of  this  modeling 
formulation  was  completed.  This  included  the  development  of  the  computational  sub-model  and 
the  model  statistics  processing  and  database  generation  software.  This  formulation  was 
implemented  within  the  production  version  of  the  CRAFT  CFD®  flow  solver  that  has  been 
developed  for  application  to  blast  and  biological  warfare  (BW)  neutralization  predictions. 
Testing  and  application  of  this  formulation  will  be  described  in  a  following  sub-section. 

6.2  Task  2.3:  Turbulence  Model  Extensions  for  Dispersed  Phase  Combustion 

Under  this  task,  a  unified  modeling  formulation  was  developed  and  implemented  to  fully 
integrate  the  modeling  formulation  developed  under  the  previous  work  period’s  efforts  with 
additional  models  to  address  particulate  combustion  and  transport  modeling.  This  effort 
consisted  of  integrating:  1)  the  spore  particle  -  turbulence  interaction  model  of  Task  1.4,  2)  a 
particle  -  turbulence  interaction  model  for  particulate  combustion,  and  3)  a  hybrid  RANS/LES 
subgrid  model.  This  development  effort  was  in  parallel  to  work  at  Georgia  Tech  to  address 
neutralization  strategies  within  the  context  of  aluminized  afterburning  munitions.  However, 
efforts  under  this  task  focused  on  large  scale,  production  level  applications  that  may  be  applied 
for  routine  analysis  needs.  With  this  unified  formulation,  the  effects  of  turbulent  flow 
fluctuations  on  the  spore  neutralization  and  particulate  combustion  may  be  captured.  Also,  this 
modeling  formulation  was  extended  to  include  a  compressible  flow  hybrid  RANS/LES  modeling 
formulation  for  turbulent  momentum  and  heat  transport.  This  hybrid  formulation  allows  for  a 
tractable,  but  advanced,  treatment  of  turbulent  transport  phenomena  within  a  large  scale 
computation.  Within  this  approach,  an  unsteady  RANS  modeling  formulation  is  employed 
within  regions  of  the  flow  where  turbulent  flow  structures  may  not  be  resolved  in  a  tractable 
manner.  These  regions  include,  for  example,  boundary  layer  flows  along  the  walls  of  a  multi¬ 
room  blast  calculations.  Turbulent  fluctuations  within  these  regions  may  not  be  tractably 
resolved  and  may  be  accurately  modeled  with  an  unsteady  RANS  formulation.  However,  away 
from  these  fine  scale  regions,  the  flow  model  transitions  to  LES  in  regions  where  large  scale 
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turbulent  structure  may  be  resolved.  Resolving  these  large  scale  structures  in  regions  of  shear 
and  mixing  will  provide  for  a  more  accurate  representation  of  the  flowfield.  Under  this  task  an 
advanced  hybrid  RANS/LES  formulation  was  coupled  with  the  spore  and  particulate  combustion 
models.  This  hybrid  model  has  also  been  extended  for  application  to  highly  compressible  flows 
that  are  characteristic  of  BW  neutralization  applications. 

The  following  sub-section  first  presents  a  brief  description  of  the  spore,  particulate  combustion, 
and  hybrid  RANS/LES  models.  This  is  followed  by  an  illustration  of  the  behavior  of  the  spore 
model  for  a  compressible  shear  layer  case,  and  a  description  of  the  preliminary  application  of  the 
modeling  approach  to  a  BW  neutralization  case. 

6.2.1  Spore  Particle  -  Turbulence  Interaction  Modeling 

As  identified  under  Task  1.2  from  the  Year  1  work  efforts,  correlations  for  the  rate  of  spore 
neutralization  may  be  specified  in  the  general  form, 

w,  =  -pY,^aT'”exp(-^/T)  (24) 

where  p  is  the  gas  phase  density,  is  the  live  spore  mass  concentration,  T  is  the  gas  phase 
temperature,  and  m,  a ,  6 ,  and  (3  are  coefficients  specified  from  experimental  correlations.  In 
this  equation,  the  spore  temperature  is  assumed  to  be  in  equilibrium  with  the  gas  phase  so  that 
the  gas  phase  temperature  appears  in  Eqn.  (24).  Application  of  Eqn.  (24)  to  turbulent  flow 
requires  characterizing  the  effects  of  spore  concentration  and  gas  phase  property  fluctuations  on 
the  mean  neutralization  rate  which  may  be  written  as, 

^  /  T )P(p,  Y^ ,  T )dpdY^dT  (25) 

where  P  is  the  joint  probability  density  function  of  the  gas  phase  density  and  temperature,  and 
the  live  spore  mass  concentration.  Given  P,  Eqn.  (25)  may  be  evaluated  to  characterize  the 
effect  of  turbulent  fluctuations  on  the  mean  neutralization  rate. 

The  pdf  in  Eqn.  (25)  may  be  specified  by  either  direct  evaluation  or  a  parameterization  strategy. 
The  direct  evaluation  approach  seeks  to  predict  the  shape  of  the  pdf  through  the  solution  of  a 
transport  equation  for  P  based  on  fundamental  principles.  This  approach  is  comprehensive,  but 
intractable  for  all  but  a  limited  range  of  simplified  problems.  Alternatively,  the  effect  of 
turbulent  fluctuations  within  Eqn.  (25)  may  be  characterized  through  a  parameterization  of  P  in 
terms  of  a  reduced  set  of  variables.  Through  such  an  approach,  a  tractable  modeling  formulation 
may  be  developed  that  may  be  applied  to  large  scale  problems  of  interest.  This  type  of 
parameterization  approach  was  chosen  to  characterize  the  effect  of  turbulent  fluctuations  on  the 
spore  neutralization  rate. 

A  parameterized  modeling  approach  has  been  developed  to  model  Eqn.  (25)  so  that  the  pdf  P  is 
parameterized  using  an  assumed  pdf  formulation  [19]  assuming  statistical  independence  between 
the  density,  temperature  and  spore  concentration.  With  this  assumption,  the  pdf  in  Eqn.  (25) 
may  be  written  as. 


P(p,TJO  «  P^(p  -  p)Pr(T-  <  T  >)Py(Y,- <  Y,  >) 


(26) 
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where  P^,  Pj,  and  Py  are  the  marginal  pdfs  the  density,  temperature  and  spore  concentration. 

In  this  equation,  the  overbar  represents  Reynolds  averaging  or  fdtering,  while  the  brackets 
represent  Favre  averaging  or  fdtering.  These  marginal  pdfs  are  then  specified  as  follows.  For 
density,  a  delta  function  is  specified  for  P^  [19].  For  and  Py,  a  beta  function  [19]  is  used 

which  is  parameterized  in  terms  of  it  first  two  moments.  For  example,  using  a  beta  function,  Py 
may  be  expressed  as. 


with. 


PyiYL)  = 


r(A)r(A) 


r(A  +  A) 


(27) 


I3.=<Y,> 


<7,>(l-<7,>) 


-1 


y^,=(l-<7,>) 


<7,>(l-<7,>) 


-1 


(28) 


where  F (x)  is  the  gamma  function.  As  seen  in  these  equations,  the  beta  function  for  7^  has  been 


parameterized  in  terms  of  the  mean  spore  concentration  (<  7^  >)  and  its  variance  (<  >).  A 

similar  expression  may  be  developed  for  Py  that  is  parameterized  in  terms  of  the  mean 
temperature  and  temperature  variance. 


With  this  parameterization  of  Eqn.  (26),  the  mean  spore  neutralization  rate  in  Eqn.  (25)  may  be 
evaluated.  With  this  formulation,  a  pdf  table  generator  code  was  developed  to  generate  a 
database  for  the  mean  spore  rate  as  a  function  of  the  mean  temperature,  temperature  variance, 
mean  spore  concentration,  and  spore  concentration  variance.  This  database  may  then  be  used  to 
specify  the  mean  spore  rate  (Eqn.  (25))  in  the  spore  transport  equation. 


This  parameterized  formulation  additionally  requires  the  specification  of  the  temperature  and 
spore  concentration  variances  that  appears  in  the  beta  representations  of  Py  and  Py .  These 

quantities  are  specified  from  addition  transport  equations.  For  the  temperature  variance,  a 
variable  turbulent  Prandtl  number  model  [20]  that  is  applicable  to  reacting  flows  and  a  hybrid 
RANS/LES  formulation  is  used.  For  the  spore  fluctuations,  a  spore  concentration  variance 
equation  was  developed  that  is  similar  to  the  species  concentration  variance  methodology  of 
Gaffney,  et  al.  [21].  This  equation  is  given  by. 


dp  <q>  <Uj><q> 


dt 


+  2.pDy 


dx, 


d 

dx, 


P(D  +  Dy) 


d  <q> 

dx. 


(29) 


5<7,> 


-  2p  ^  ^  ^  -2<Yy>(l)y+  2Y,(l), 


L^L 


J  / 


where  <  q  >=<  7^  > ,  and  ty  is  the  turbulent  mixing  time  scale  as  specified  from  the  hybrid 


RANS/LES  model.  The  term  YyCOy  is  also  modeled  using  the  assumed  pdf  methodology  of  Eqn. 

(26).  The  parameterized  model  of  this  term  is  also  stored  within  the  model  database  that  is 
generated  by  the  pdf  table  generator. 
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6.2.2  Particle  -  Turbulence  Interaction  Modeling  for  Aluminum  Particulate 
Combustion 

For  turbulent  flow,  velocity  and  composition  fluctuations  will  influence  the  burning  rate  of  solid 
aluminum  particulates  within  afterburning  munitions  that  may  be  used  for  spore  neutralization. 
To  account  for  these  fluctuations,  two  approaches  may  be  considered.  The  first  is  a  microscopic 
approach  that  seeks  to  model  the  direct  interaction  of  the  flow  with  the  flame  structure 
surrounding  the  particles  from  a  first  principles  perspective.  The  second  is  a  macroscopic 
approach  that  seeks  to  account  for  flowfield  fluctuations  on  the  modeled  laminar  burning  rate. 
The  microscopic  approach  is  more  fundamental,  but  intractable  from  the  perspective  of  applying 
it  routinely  to  problems  of  practical  interest.  The  macroscopic  approach  is  more  pragmatic  and 
empirical,  but  provides  a  viable  alternative  to  develop  a  tractable  formulation.  As  a 
consequence,  modeling  efforts  focused  on  the  development  of  a  macroscopic  particle  - 
turbulence  interaction  model  that  may  be  routinely  applied  to  problems  of  interest. 

For  the  macroscopic  approach,  the  model  sought  to  account  for  the  effects  of  turbulent 
fluctuations  on  the  modeled  laminar  particle  burning  rate.  Within  CRAFT  Tech’s  Eulerian 
dispersed  phase  particle  combustion  modeling  formulation,  the  laminar  aluminum  particle  mass 
density  consumption  rate  is  expressed  as. 


mp  =  -IJirlNp^  (30) 

at 

where  Pp  is  the  particle  mass  density,  pm  is  the  material  density,  N  is  the  particle  number  density, 
and  rp  and  dp  are  the  particle  radius  and  diameter,  respectively.  The  time  derivative  in  Eqn.  (30) 
is  also  expressed  as. 


^dp _ 1 _ 

dt  ndl'^Ap  +  mdp-^A^ 

with, 


(31) 


A=- 


1 


1 


4  exp(-E^  /  RT  )X..  1  -  £ 


A.  = 


1 


ro-2po.i^^^  i-f" 


(32) 

(33) 


where  n  =  0.3,  m  =  1.8,  A]  =  5.5x10^  A2  =  7.35x10'*’,  e  =  0.05,  Eb  =  73.6  kJ/mol,  R  is  the  gas 
constant,  Tp  is  the  particle  surface  temperature,  P  is  the  pressure,  T  is  the  gas  temperature.  The 
variable  Xef/is  the  sum  of  gas  phase  mole  fractions  of  O2,  H2O  and  CO2  as  given  by, 

Xj  -  +  0.22X0.  (34) 


Eqns.  (30)  -  (34)  were  developed  from  the  Yetter’s  hybrid  combustion  model  [22]  for  aluminum 
particle  burning.  This  model  formulation  accounts  for  both  the  diffusion-controlled  and 
kinetically-controlled  burning  limits,  and  a  transition  between  the  two.  With  this  formulation, 
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the  particle  mass  density  consumption  rate  in  Eqn.  (30)  may  be  written  as  a  function  of  the 
particle  and  gas  phase  properties  as, 


m,-f(r^.T,.P.T.XJ  (35) 

where  Xu  are  the  species  mole  fractions.  This  equation  is  a  representation  of  the  laminar  burning 
rate  model  within  CRAFT  Tech’s  dispersed  phase  combustion  model  for  aluminum  particles. 

With  the  laminar  rate  model  given  by  Eqn.  (35),  a  macroscopic  level  model  for  the  effect  of 
turbulent  fluctuations  on  the  rate  expression  could  be  represented  by, 

^ -Jpdf(r^.T^.P.T,X,)f(r^.T^,P,T,XJdr/T^dPdTdX,  (36) 

where  pdf  ( rp,Tp,P,T ,X^)  is  the  joint  probability  density  function  of  the  particle  and  gas  phase 

properties  and  the  overbar  represents  Reynolds  averaging.  With  Eqn.  (36),  a  parameterized 
model  for  this  mean  rate  expression  has  been  developed.  As  discussed  earlier,  a  parameterized 
modeling  formulation  seeks  to  represent  the  model  statistics  in  terms  of  a  reduced  set  of 
variables.  These  model  statistics  may  then  be  stored  within  a  database.  A  flow  solver  may  then 
retrieve  the  model  statistics  from  the  database  in  a  fast  and  efficient  manner. 

To  develop  a  parameterized  model  for  Eqn.  (36),  it  should  first  be  recognized  that  the  particle 
and  gas  phase  properties  for  dispersed  phase  combustion  are  effectively  uncorrelated.  That  is, 
specific  gas  phase  properties  do  not  necessary  correspond  to  specific  particle  phase  properties. 
For  example,  high  instantaneous  values  of  gas  phase  temperature,  T,  do  not  necessarily 
correspond  to  high  instantaneous  values  of  particle  temperature,  or  specific  values  of  particle 
radius.  The  particle  and  gas  phase  properties  will  only  be  strongly  correlated  when  the  two 
phases  are  in  equilibrium.  As  a  result,  the  particle  and  gas  phase  properties  may  be  assumed  to  be 
statistically  independent.  With  this  assumption,  the  joint  pdf  in  Eqn.  (36)  may  be  expressed  as  a 
combination  of  the  marginal  pdfs  of  these  variables  as, 

pdf(r^.T^.P.T.X,).pdf(r^.T;>pdf(P.T.Xp>  (37) 

With  this  representation  of  the  joint  pdf,  an  assumed  pdf  model  [19]  may  be  used  to  develop  an 
approximation  for  the  RHS  of  Eqn.  (37).  To  develop  a  first  order  model,  the  pdf  of  the  particle 
properties  is  assumed  to  be  composed  of  delta  functions  at  the  mean  particle  values.  For  the  gas 
phase  pdf,  the  pressure  is  also  assumed  to  be  statistically  independent  of  the  temperature  and 
species.  These  assumptions  are  reasonable  for  the  present  application  because  the  particle  and 
gas  phase  properties  are  effectively  uncorrelated  as  discuss  earlier.  With  these  assumptions, 
Eqn.  (37)  may  be  rewritten  as. 


pdf(r^  .T^.P.T.X,).dO-^-r^)d(T^-T^  )pdf(P)pdf(T,  XJ  (38) 

For  the  gas  phase  pdfs,  assumed  pdf  methods  typically  employ  a  delta  function  for  the  marginal 
pdf  of  pressure.  The  task  then  becomes  to  develop  a  parameterized  model  for  pdf  (T ,X,^) .  For 
afterburning  munitions,  particle  -  turbulence  interactions  will  primarily  occur  in  the  shear  layer 
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region  between  the  ambient  fluid  and  the  blast  constituents  where  mixing  and  combustion  may 
occur.  Under  these  conditions,  the  mixture  gas  phase  properties  could  be  parameterized  in  terms 
of  mixture  fraction  (Z)  and  scalar  dissipation  (x),  similar  to  the  parameterization  used  by 
Sankaran,  et  al.  [23]  for  gas  phase  combustion  and  for  flamelet  models  [19].  Consequently,  the 
local  temperature  and  species  within  the  shear  layer  may  be  written  as  T  =  T(Z,x,P)  and 

=  Xj^(Z,x,P),  with  the  pressure  being  included  due  to  compressibility  effects  within  the 

flow.  This  parameterization  of  the  gas  phase  properties  may  be  accomplished  using  a  flamelet 
model  [19]  given  representative  thermodynamic  properties  of  the  flow  (i.e.,  fuel  and  oxidizer 
composition,  and  pressure  range).  With  this  representation,  and  the  aforementioned  assumption 
regarding  the  statistical  independence  of  the  pressure  and  temperature  and  species,  Eqn.  (38) 
then  becomes, 

pdf(r^.T^.P. T.X,)~ S(r^  - 7^)S(T^  -  %)S(P - P)&(x-  <  X  >)pdf(Z)  (39) 

where  the  typical  assumptions  regarding  scalar  dissipation  from  assumed  pdf  methods  have  been 
employed  [19]. 

Within  Eqn.  (39)  ,  the  pdf  of  mixture  fraction  may  be  approximated  using  a  beta  function  [19] 
The  mean  particle  density  consumption  rate  in  Eqn.  (39)  may  be  evaluated  as  a  function  of  the 
mean  flow  and  particle  properties  as, 

^  =  g(y,,f'^)h(P,<  X  >,<  Z  >,<  Z'^  >)  (40) 

where  the  functions  g  and  h  are  independent  due  to  the  previous  assumption  that  the  particle  and 
flow  properties  are  uncorrelated.  Alternatively,  a  function  F  may  be  defined  as  the  ratio  of 
turbulent-to-laminar  burning  rate  as. 


F(r„,T„,P,<  X>,<Z>,<Z'^  >,<T  >)  = 


m 


(41) 


'p,lam 


where  is  the  consumption  rate  defined  by  Eqn.  (35)  evaluated  based  on  mean  flow 

quantities  only.  The  term  is  the  value  of  the  consumption  rate  that  neglects  all  turbulent 

fluctuations,  or  that  is  the  effective  laminar  rate. 


With  Eqns.  (39)-(41),  a  database  for  the  function  F  may  be  evaluated  given  the  representative 
thermodynamics  conditions  of  the  problem.  This  database  may  then  be  deployed  within  a  CFD 
flow  solver  so  that  the  mean  particle  density  consumption  rate  may  be  evaluated  as. 


•  m 


p,lam 


F( r„,T„,P,<  X  >-<  >.<  T  >) 


(42) 


By  defining  the  mean  consumption  rate  in  this  manner,  the  first  order  effect  of  turbulent 
flowfield  fluctuations  may  be  included.  To  apply  this  model,  the  CFD  flow  solver  must 
additionally  solve  standard  transport  equations  for  <  Z  >  and  <  Z'^  > .  These  transport  equation 
may  be  specified  using  the  compressible  flow,  variable  turbulent  Schmidt  number  formulation  of 
Brinckman,  et  al.  [24]. 
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With  this  modeling  formulation,  the  CFD  flow  solver  is  only  required  to  include  just  two 
additional  transport  equations.  The  multiplying  function  F  is  also  retrieved  from  a  database  in  a 
fast  and  efficient  manner.  As  a  result,  this  first  order  turbulence  -  particle  interaction  model  is 
computationally  inexpensive  to  employ  for  large  scale  applications. 


6.2.3  Hybrid  RANS/LES  Model  for  Momentum  and  Heat  Transport 

The  hybrid  RANS/LES  (HRLES)  model  employs  the  standard  k-8  equations  of  a  typical  RANS 
formulation.  However,  the  eddy  viscosity  is  suitably  scaled  down  in  EES  regions  based  on  an 
assessment  of  the  local  resolution  levels  and  the  local  range  of  scales  that  this  resolution  would 
permit.  This  assessment  is  based  on  the  RANS  estimate  for  the  turbulent  kinetic  energy,  and  an 
estimate  of  the  range  of  turbulent  scales  that  could  be  supported  on  the  computational  grid, 
assuming  a  pre-specified  empirical  form  of  the  turbulent  kinetic  energy  (TKE)  spectrum.  This 
formulation  is  based  on  the  model  of  Arunajatesan  and  Sinha  [25],  which  was  developed 
assuming  low  speed  flow.  In  this  regard,  the  Karman-Pao  energy  spectrum  model,  which  is 
widely  used  and  well  accepted  as  a  good  candidate  for  purely  incompressible  flows,  was  used  to 
develop  the  original  formulation  [25].  However,  this  formulation  has  been  upgraded  to  account 
for  compressibility  effects  to  extend  its  range  of  application  to  blast  predictions. 


Flow  ccompressibility  affects  the  TKE  spectrum  as  illustrated  in  Figure  3  [26].  From  this  figure, 
note  that  both  the  slope  of  the  spectrum  in  the  inertial  wavenumber  regime,  as  well  as  the  decay 
rate  in  the  dissipation  range  are  affected  by  compressibility.  The  non-dimensional 
compressibility  parameter  (used  to  characterize  the  magnitude  of  compressibility)  that  is  widely 

used  in  the  literature  is  that  of  the  turbulent  Mach  number,  ,  which  is  a  ratio  of  the 

^  \yRT 


turbulent  velocity  scale  to  the  acoustic  speed. 


To  account  for  the  effects  of  compressibility,  the  turbulent  kinetic  energy  is  decomposed  into  a 
solenoidal  and  a  compressible  component,  written  as: 


(43) 


where. 


4  r 


1  + 


f-1 


-17/ 

2  1  /6 


exp 


/  3  4/\ 

2 


(44) 


is  the  solenoidal  component  that  is  represented  by  the  Karman-Pao  empirical  model,  and. 


(45) 
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is  the  compressible  part  of  the  energy  that  scales  with  the  square  of  turbulent  Mach  number,  . 
Note  that  these  equations  are  non-dimensionalized  and  use  the  non-dimensional  wavenumber, 
k=kri,  where  ri  =  \y  le\  is  the  Kolmogorov  dissipation  length  scale.  The  function 

parameterizing  the  compressible  component  of  the  spectrum  is  curve-fit  from  data  (as  seen  in 
Figure  4)  in  terms  of  and  k^ .  The  ratio  of  the  compressible  energy  to  the  solenoidal  energy 
(see  Figure  5)  is  also  curve-fit  as  a  function  of  and  the  Taylor-scale  Reynolds  number,  given 
by 


Re^  = 


(46) 


Figure  3.  Turbulent  Kinetic  Energy  Spectra 
[26] 


Figure  4.  Compressible  Energy  Spectra 
Normalized  by  [26] 
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Figure  5.  Acoustic  Energy  Divided  by  the  Solenoidal  Energy  as  a  Function  of  [26] 


Note  that  the  model  parameters  Cs,  Cc  and  need  to  be  determined  in  terms  of  the  two 
independent  variables  Rex  and  M^.  Cs  and  Cc  are  the  constants  multiplying  the  energy 
components  (solenoidal)  and  (compressible)  respectively,  while  is  the  dominant 

energy  carrying  wavenumber  (non-dimensional)  in  the  spectrum.  Just  as  in  the  original  model 
formulation  [25],  the  eddy-viscosity  in  the  LES  regions  is  damped  by  hybrid  function  that  is 
now  a  function  of  both  Rex  and  Mr.. 


V. 


T,LES 


f hybrid^'. 


hybrid  T,RANS 


_  r  f  rans 

J  hybrid  J  /li  fi 


(47) 


In  order  to  determine  the  three  model  parameters  G,  Cc  and  k^ ,  the  following  set  of  equations 
(three)  need  to  be  solved: 


Re,.^y[C,£,h)  +  C,M/£,(i) 


(48) 


The  first  equation  above  is  that  the  sum  of  the  integrals  of  the  two  components  of  the  energy 
(non-dimensional)  which  should  equal  Rex,  while  the  second  equation  is  that  the  sum  of  the 
integrals  of  two  components  of  the  dissipation  of  energy  that  is  equal  to  the  non-dimensional 
dissipation  rate,  e  =  1.  The  last  equation  is  a  ratio  of  the  compressible  to  the  solenoidal  energy. 
The  above  set  of  three  equations  is  solved  for  the  three  model  parameters  for  a  range  of  values  of 


Final  Report 


27 


CRAFT  Tech  Report  C496 


the  independent  variables  Rex  and  M^.  The  surfaces  of  these  model  parameters  are  shown  in 
Figure  6.  Note  that  all  the  parameters  asymptote  to  a  constant  value  in  the  high  Reynolds 
number  limit.  Once  the  three  model  parameters  are  determined,  the  entire  energy  spectrum  is 
known.  Then,  for  different  values  of  the  mesh-scale  to  the  energy-carrying  length  scale,  A  /  4  the 

subgrid-scale  eddy-viscosity,v^^^^is  evaluated  by  integrating  the  energy  and  the  dissipation 
spectrum  from  the  mesh-length  scale.  A,  to  infinity  (or  the  smallest  scales,  i.e.,  the  Kolmogorov 
scale,  rj  =  \v  ! e\  ). 


=f^^E{k)dk-  =f^^D(k)dk 
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The  surface  of  the /hybrid  function  is  also  shown  in  Figure  6.  These  surface  functions  are  curve-fit 
and  explicitly  defined  in  terms  of  the  two  independent  variables  Rex  and  Mr..  Hence,  at  any  given 
point  in  the  flow  field,  the  determination  of  these  quantities  along  with  a  definition  of  the  local 
mesh-length  scale  can  be  used  to  evaluate  the  /hybrid  function.  This  essentially  damps  the  eddy- 
viscosity  contribution  for  the  given  mesh-length  scale  based  on  that  portion  of  energy  un¬ 
sustainable  by  the  local  mesh,  and  hence  is  the  “unresolved”  energy  local  to  this  point. 

The  compressibility  corrected  (c.c.)  HRLES  model  was  first  tested  on  a  hot,  over-expanded 
Mach  1.5  jet.  Figure  7  presents  instantaneous  Mach  number  plots  comparing  the  incompressible 
HRLES  model  and  the  HRLES  model  with  c.c.  A  limitation  of  the  incompressible  model  was  a 
general  over-prediction  of  eddy  viscosity  in  the  initial  shear  layer  of  the  jet.  This  eddy  viscosity 
was  calculated  assuming  the  incompressible  energy  spectra  and  therefore  the  original  model 
tended  to  smooth  the  shear  layer  and  delayed  the  evolution  of  jet  turbulence.  As  can  be  seen  in 
the  lower  part  of  this  figure,  correcting  the  energy  spectra  for  compressibility  naturally  reduces 
the  amount  of  eddy  viscosity  and  allows  the  jet  shear  layer  to  “trip”  sooner  (closer  to  the  nozzle 
lip).  The  amount  of  eddy  viscosity  reduction  can  clearly  be  seen  in  Figure  8. 
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Figure  6.  Surfaces  of  Model  Parameters  C*,  Cc,  and  i\itfhybrid  Damping  Function  for  the 

Eddy-viscosity  in  the  HRLES  c.c  Model 


Figure  7.  Mach  Number  Comparison  of  Incompressible  HRLES  Model  with 
Compressibility  Corrected  HRLES  Model 
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Figure  8.  Eddy  Viscosity  Ratio  Comparison  of  Incompressible  HRLES  Model  with  c.c. 

HRLES  Model 

6.2.4  Model  Testing  and  Evaluation 

As  discussed  earlier,  this  unified  formulation  was  implemented  within  the  production  version  of 
the  CRAFT  CFD®  flow  solver  that  has  been  developed  for  applieation  to  blast  and  biologieal 
warfare  (BW)  neutralization  predietions.  As  an  initial  test  of  the  formulation,  a  high  speed  shear 
layer  eonfiguration  was  eonsidered  to  evaluate  the  effects  of  the  turbulenee  modeling  on  live 
spore  depletion  or  neutralization.  This  configuration  was  chosen  for  this  initial  test  beeause 
spore  neutralization  within  a  blast  configuration  will  primarily  occur  in  turbulent  shear  regions 
between  the  ambient  conditions  (containing  the  spore  eonstituents)  and  the  blast  flow.  For  this 
shear  layer  configuration,  the  low  speed  flow  side  eonsisted  of  air  at  1  atm  of  pressure  and  an 
initial  temperature  of  300  K,  and  with  a  flow  Maeh  number  of  1.2.  This  low  speed  flow  was  also 
contaminated  with  spore  eonstituents  at  a  eoncentration  of  10  ppm.  The  properties  of  these 
hypothetieal  spores  with  respeet  to  the  neutralization  rate  in  Eqn.  (24)  were  ctr  =  100,  6  =  \,m  = 
0,  and  =  350K.  For  the  high  speed  side,  the  blast  flow  was  modeled  as  reaction  products  with 
an  excess  of  CO  and  H2  in  equal  proportions  at  a  pressure  of  1  am  and  temperature  of  2000  K, 
and  with  a  flow  Maeh  number  of  2.5. 

With  these  conditions,  the  spore  model  pre-processor  developed  under  this  program  was  used  to 
generate  the  database  for  the  spore  particle  -  turbulent  interaetion  model.  This  database  includes 
model  data  for  the  mean  or  filtered  live  spore  destruction  rate,  co^ ,  and  the  source  term  in 

the  spore  eoncentration  varianee  equation  given  in  Eqn.  (29).  These  two  terms  are  parameterized 
within  the  database  as. 


=  f  (<¥,>,<  Y'l>,<T>,<T'^>)  (50) 

Figure  9  illustrates  the  effeet  of  the  pdf  distribution  used  to  model  the  temperature  fluctuations 
on  the  mean  spore  destruetion  rate.  This  figure  presents  contours  of  the  ratio  of  the  turbulent 
spore  destruetion  rate  to  the  laminar  destruetion  rate  that  neglects  temperature  fluetuations. 
These  eontours  are  plotted  in  the  mean  temperature  -  temperature  variance  plane.  In  this  figure, 


Final  Report 


30 


CRAFT  Tech  Report  C496 


the  upper  boundary  in  temperature  variance  is  nonlinear  due  to  realizability  constraints  on  the 
assumed  pdf  distribution  for  the  temperature.  As  seen  in  the  figure,  the  maximum  realizable 
temperature  variance  is  equal  to  zero  at  the  minimum  and  maximum  specified  temperature  range, 
and  reaches  a  maximum  value  in  between  these  limits  at  ~  500  K.  Over  this  realizable  range,  the 
effect  of  temperature  fluctuations  as  seen  in  this  figure  is  to  reduce  the  mean  spore  destruction 
rate.  As  the  temperature  fluctuations  become  large,  the  mean  destruction  rate  is  seen  to  be 
reduced  by  ~  50%.  This  model  prediction  has  important  implications  for  the  assessment  of 
prospective  neutralization  strategies.  As  the  results  in  Figure  9  indicate,  neglecting  the  effect  of 
temperature  fluctuations  on  the  destruction  rate  could  significantly  over  estimate  the 
effectiveness  of  a  neutralization  strategy.  As  a  result,  it  is  important  to  include  the  effect  of  these 
fluctuations  on  any  predictive  assessment  to  provide  a  more  conservative  estimate  of 
neutralization  effectiveness. 


Figure  9.  Contours  of  the  ratio  of  the  turbulent-to-laminar  live  spore  destruction  rate  as  a 
function  of  mean  temperature  and  normalized  temperature  variance 

Results  for  the  application  of  the  spore  interaction  model  to  the  high  speed  shear  layer  case  are 
presented  in  Figure  10.  This  figure  presents  contours  of  the  normalized  spore  concentration 
within  the  shear  layer  for  the  laminar  rate  model  and  the  HRLES  -  spore  interaction  model.  In 
this  figure,  the  low  speed  side  of  the  shear  layer  with  the  spore  contamination  is  on  the  top  half, 
while  the  high  speed  blast  flow  is  on  the  bottom  half  For  this  case,  the  temperature  variance  for 
the  spore  interaction  model  case  has  been  set  to  the  maximum  realizable  value  to  illustrate  the 
maximum  possible  effect  of  the  model.  Also,  the  streamwise  scale  of  these  contour  plots  has 
been  dramatically  scaled  for  clarity.  From  Figure  10,  it  is  clear  that  under  these  conditions  that 
spore  destruction  within  the  shear  layer  has  been  substantially  reduced.  The  spore  concentration 
in  this  figure  is  also  seen  to  have  not  reached  self-similarity  within  the  shear  layer,  as  indicated 
by  the  nonlinear  nature  of  the  contours.  As  a  result,  the  laminar  and  turbulent  rate  model  results 
will  continue  to  diverge  further  downstream.  The  lack  of  self-similarity  in  this  figure  is  a  result 
of  that  fact  that  the  chemical  reactions  within  the  shear  layer  have  not  reached  self-similarity 
over  the  spatial  domain  of  this  figure.  For  this  case,  excess  fuel  in  the  blast  flow  is  igniting  along 
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the  shear  layer  and  the  ehemical  reaetions  have  not  yet  reaehed  a  steady  state.  Consequently,  the 
results  presented  in  Figure  10  also  indicate  that  the  spore  interaction  model  is  also  sensitive  to 
ignition  phenomena  occurring  within  the  gas  phase. 


With  this  test  of  the  spore  interaction  model  completed,  the  model  implementation  within 
CRAFT  CFD®  was  considered  verified. 


(a)  Laminar  spore  destruction  rate  model  (b)  HRLES  -  spore  interaction  model 

Figure  10.  Contours  of  the  normalized  spore  concentration  within  the  reacting  high  speed 

shear  layer 


6.2.5  Biological  AD  Application 

With  the  spore  interaction  model  implementation  now  verified,  the  formulation  was  then  applied 
to  a  high-fidelity  blast  event  simulation  of  a  representative  aluminized  high  explosive  (HE) 
charge  in  a  generic  configuration  featuring  an  inner  room  within  an  outer  room  communicating 
through  various  openings.  For  this  demonstration,  the  biological  agent  consisting  of  b.T.  spores 
was  assumed  to  be  contained  within  8  canisters  distributed  in  a  semi-circular  pattern  centered 
around  the  HE  charge  in  the  inner  room.  Due  to  the  significant  computational  cost  of  time- 
accurate  3-D  high-fidelity  simulations  of  blast/ AD  scenarios,  the  biological  agent  release  was 
intentionally  localized  close  to  the  floor  of  the  inner  room,  where  turbulent  fluctuations  in 
temperature  and  the  establishment  of  a  high- temperature  gas  environment  (HTGE)  in  direct 
contact  with  the  spores  are  most  likely  to  occur.  This  initialization  is  aimed  at  enhancing  the 
effect  of  the  turbulent  HTGE  on  the  spores  in  the  very  early  stages  of  the  post-detonation  and 
therefore  at  reducing  the  required  duration  of  the  simulated  overall  event. 

When  the  HE  charge  is  first  detonated,  the  resulting  blast  front  propagates  outwards  in  a  complex 
three-dimensional  pattern.  Shocks  are  reflected  off  the  floor  and  later  from  the  walls.  During  this 
process,  blast  afterburning  takes  place  as  the  detonation  products  mix  with  the  air  available  in  the 
inner  room.  As  a  result,  the  heat  release  from  the  blast  afterburning  further  raises  the  temperature 
inside  the  inner  room,  thus  providing  a  suitable  HTGE  for  the  spore  demise.  Due  to  the  blast 
over-pressure,  all  vents  and  the  door  are  at  first  choked  and  only  outflow  from  the  inner  room 
takes  place.  As  the  simulation  progresses,  it  is  expected  that,  after  reaching  the  peak  pressure 
level  in  the  inner  room,  the  pressure  starts  dropping.  As  blast  afterburning  takes  place,  turbulent 
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mixing  in  a  complex  highly  three-dimensional  pattern  eontrols  and  limits  the  rate  of  combustion. 
When  pressure  equilibration  between  the  two  rooms  is  reaehed,  the  vents  and  the  doorway  are  no 
longer  choked  and  flow  reversal  will  take  place.  Fresh  air  from  the  outer  room  enters  the  inner 
room  where  reaetants  are  still  available  at  oxygen-lean  eonditions  and  hot  products  exit  from  the 
inner  room  in  a  pulsating  flow  pattern. 

The  b.  T.  spores  are  introduced  into  the  inner  room  of  the  geometry  during  the  HE  charge  post¬ 
detonation  with  a  delay  representative  of  the  time  required  by  HE  eharge  casing  fragments  to  hit 
and  rupture  the  spore  eontainers  and  release  the  elouds  of  biological  agent.  Based  on 
representative  empirical  correlation  for  spores,  the  temperature  dependenee  of  the  turbulent  spore 
demise  rate  is  illustrated  in  Figure  11.  Speeifically,  Figure  11(a)  presents  a  contour  plot  of  the 
turbulent  rate  in  temperature-temperature  variance  space,  indicating  a  strong  decay  in  the  rate 
with  decreasing  temperature  irrespeetive  of  the  level  of  turbulent  fluetuations.  On  the  other  hand. 
Figure  11(b)  and  Figure  11(c)  show  the  ratio  of  laminar  to  turbulent  spore  demise  rate  on  two 
different  seales.  This  is  intended  to  highlight  the  moderately  inhibiting  effect  of  turbulent 
fluetuations  on  the  rate  at  high  temperature  (see  Figure  1 1(b))  and  the  major  enhaneing  effect  in 
the  lower  temperature  range.  However,  since  the  rate  is  small  at  low  temperature,  the  overall 
effeetiveness  of  the  rate  enhaneement  is  mitigated.  What  these  plots  indieate  is  that  a  visible 
inerease  in  the  spore  demise  rate  ean  be  attained  in  the  mid  temperature  range,  if  signifieant 
turbulenee  levels  are  present. 


Temperature  Temperature  Temperature 


(a)  (b)  (e) 

Figure  11.  Contour  plots  in  temperature  -  temperature  variance  space  of  (a)  turbulent 
spore  demise  rate,  and  (b-c)  ratio  of  laminar  to  turbulent  spore  demise  rates  using  two 
different  scales  to  emphasize  enhancing/inhibiting  effects. 

A  HRLES  simulation  on  a  moderately  coarse  grid  based  on  half-plane  symmetry  and  activating 
the  spore  partiele-turbulence  interaetion  model  was  carried  out.  Due  to  the  signifieant 
computational  cost,  numerical  stiffness  and  restrietive  integration  time  steps,  the  early  post¬ 
detonation  dynamics  of  fireball  expansion  was  simulated  for  a  duration  of  8  ms.  The  evolution  in 
time  of  the  volume-integrated  normalized  mass  of  alive  spores  is  presented  in  Figure  12, 
indieating  a  rapid  initial  demise  of  the  spores  as  the  shoek  wave  goes  through  the  cloud  of 
biological  agent  and  a  reduetion  in  the  rate  as  more  uniform  HTGE  is  established  in  its 
proximity. 
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Figure  12.  Volume-integrated  normalized  alive  spore  mass  as  a  function  of  simulation  time. 


Selected  snapshots  of  the  pressure  distribution  inside  the  inner  room  are  shown  in  Figure  13.  The 
alive  spore  cloud  is  identified  by  iso-surfaces  shown  in  yellow  and  the  hot  plume  venting  into  the 
outer  room  from  the  ceiling  vents  and  doorway  by  temperature  iso-surfaces  in  orange.  Similarly 
the  temperature  distribution  corresponding  to  the  same  instants  in  time  is  presented  in  Figure  14. 
The  deformation  and  destruction  of  the  cloud  of  alive  spores  can  be  observed  as  the  shock  wave 
propagates  through  the  cloud  and  gets  reflected  off  the  inner  room  walls  and  ceiling.  Due  to  the 
way  the  HE  charge  was  detonated,  resulting  in  a  predominantly  downward-directed  blast 
propagation,  significant  HTGE  is  established  along  the  floor  of  the  room,  as  indicated  by  the 
contours  in  Figure  14.  While  this  simulation  could  potentially  be  carried  out  for  hundreds  of 
milliseconds  to  capture  the  entire  dynamics  of  the  blast/ AD  scenario,  the  present  simulation  of 
the  early  post-detonation  stages  serves  as  a  demonstration  of  the  upgraded  modeling  framework. 
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Figure  13.  Pressure  distribution  in  inner  room  for  selected  instants  in  time  with  alive  spore 
iso-surface  shown  in  yellow  and  temperature  iso-surface  shown  in  orange  (outer  room 

only). 
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Figure  14.  Temperature  distribution  in  inner  room  for  selected  instants  in  time  with  alive 

spore  iso-surface  shown  in  yellow. 
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6.3  Task  2.4:  Turbulence  Model  Extensions  for  the  Continuum  Phase  Combustion 

This  task  considers  turbulence  model  extensions  to  include  turbulence  -  chemistry  interaction 
modeling  for  the  gas  phase  combustion.  The  focus  of  this  modeling  effort  is  to  develop  and 
implement  turbulent  combustion  models  that  are  tractable  for  large  scale,  production  level  blast 
applications.  As  a  result,  the  first  phase  of  this  task  considers  the  implementation  of  a  standard, 
generally  applicable,  first  order  turbulent  combustion  model.  The  second  phase  of  this  task 
considers  extension  of  a  more  comprehensive  modeling  approach  for  blast  applications.  This 
work  period,  phase  one  of  this  task  was  completed  and  phase  two  was  initiated  as  will  be 
describe  next. 

For  the  first  phase  of  this  task,  the  Eddy  Dissipation  Concept  (EDC)  [27]  model  was 
implemented  within  CRAFT  CFD®.  This  formulation  was  originally  developed  by  Magnussen 
and  Hjertager  [27]  to  provide  a  first  order  effect  of  turbulent  fluctuations  on  chemical  reactions. 
The  EDC  model  recognizes  that  combustion  only  occurs  in  molecularly  mixed  fine  scale 
structures  of  the  flow.  These  fine  structures  only  comprise  a  small  fraction  of  the  total  volume  of 
the  fluid.  Geometrically,  these  structures  have  a  thickness  on  the  order  of  the  Kolmogorov 
microscale  and  form  broad  sheets  that  are  convoluted  by  large  scale  turbulence.  This  theory 
assumes  these  fine  scale  structures  are  thick  compared  with  the  Kolmogorov  microscale  and  are 
therefore  homogeneously  mixed  by  small  scale  turbulence.  These  structures  may  then  be 
modeled  using  a  perfectly  stirred  reactor  (PSR).  When  chemical  kinetic  rates  are  fast, 
combustion  in  the  fine  scale  structures  depends  only  on  the  mass  transfer  rate  from  the 
surrounding  inert  fluid  to  the  fine  structures.  When  chemical  kinetic  rates  are  slow,  combustion 
within  the  fine  structures  will  depend  on  the  fine  scale  Damkohler  number,  Da  *,  given  by. 


Da*=Tlt^Uen,  (51) 

where  t  is  the  fine  structure  time  scale  or  small  scale  turbulent  strain  rate.  For  RANS 
applications,  Byggstoyl  and  Magnussen  [28]  developed  the  following  expression  for  t  based  on 
turbulent  scaling  relations. 


r  =  C, 


(52) 


where  e  is  the  turbulent  kinetic  energy  dissipation  rate,  v  is  the  mean  kinematic  viscosity  and  Cs 
is  a  calibration  coefficient  typically  take  as  G  =  0.411.  The  volume  fraction  of  the  reacting  fine 
scales,  ,  is  also  given  by  [29], 
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(53) 


where  k  is  the  turbulent  kinetic  energy  and  is  a  calibration  coefficient  typically  taken  as 

Cj.  =2.13.  With  Eqn.  (52)  and  (53),  the  EDC  model  may  be  implemented  for  finite  rate 

chemistry  using  an  stiff  ordinary  differential  equation  (ode)  chemistry  solver.  In  this  context,  the 
mean  or  filtered  chemical  chemical  production  rate  for  the  kth  species  is  given  by, 
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(54) 
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The  primary  effect  of  the  EDC  model  is  to  reduce  chemical  production  when  turbulent 
fluctuations  are  present.  This  typically  has  the  effect  of  reducing  flame  temperatures  and  product 
formation.  For  example,  Figure  15  presents  radial  temperature  profdes  for  a  low  speed  hydrogen 
-  air  jet  diffusion  flame  [30]  at  an  axial  station  ofx/D  =  22.5.  Included  in  this  plot  are  the  FDC 
model  results,  laminar  chemical  rate  results,  and  experimental  data.  As  seen  in  this  figure,  the 
EDC  model  reduces  the  peak  mean  temperature  to  produce  a  slightly  better  prediction  than  when 
turbulent  fluctuations  are  neglected  within  the  chemical  production  rate  formulation. 


y/D 


Figure  15.  Favre  mean  temperature  profile  at  x/D  =  22.5. 


As  discussed  earlier,  the  EDC  model  is  generally  applicable  and  computationally  inexpensive, 
but  it  only  provides  for  a  first  order  accurate  representation  of  turbulence  -  chemistry 
interactions  within  the  flow.  For  phase  two  of  this  task,  efforts  were  initiated  this  work  period  to 
enhance  the  accuracy  of  the  turbulent  combustion  modeling  while  still  maintaining  a  tractable 
formulation.  In  coordination  with  Georgia  Tech,  the  focus  of  these  efforts  was  on  employing  the 
linear-eddy  model  (LEM)  [31]  [32]  for  turbulent  combustion.  The  LEM  is  a  comprehensive 
stochastic  mixing  model  that  separately  treats  molecular  diffusion  (with  finite-rate  kinetics)  and 
small  scale  turbulent  stirring.  This  modeling  formulation  has  been  the  focus  of  Georgia  Tech 
combustion  modeling  for  many  years  due  to  the  advanced  nature  of  the  model.  The  LEM  has 
been  implemented  by  Georgia  Tech  as  an  LES  subgrid  model  [33].  This  implementation 
includes  a  LEM  stochastic  simulation  within  each  LES  subgrid  domain  during  a  run-time 
simulation.  Though  very  comprehensive  and  accurate,  this  implementation  is  intractable  for 
large  scale,  production  level  blast  applications  due  to  the  computational  cost  of  the  inline 
stochastic  simulations.  As  an  alternative,  the  LEM  may  be  included  through  a  procedure  called 
stochastic  model  parameterization  [34].  Using  this  procedure,  statistics  from  the  LEM  are  pre- 
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computed  before  a  run-time  simulation  and  parameterized  in  terms  of  a  redueed  set  of  variables 
to  form  a  computationally  inexpensive  database  of  the  model.  This  type  of  approaeh  has  been 
developed  for  nonpremixed,  partially  premixed  [34],  and  premixed  [35]  flows  in  the  eontext  of 
the  low  speed  flow  regime.  This  type  of  approach  has  the  advantage  of  eomputational  efficieney 
and  eould  be  a  viable  modeling  approaeh  for  BW  neutralization  predictions.  However,  to  be 
feasible  this  formulation  will  require  its  extension  to  the  compressible  flow  regime. 
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7.0  YEAR  3  AND  4  ACTIVITIES 


Due  to  the  sensitive  nature  of  the  work  performed  in  Year  3  and  4  pertaining  to  applications  to 
agent  defeat  problems  of  interest  to  DTRA,  two  (2)  separate  detailed  reports  were  delivered 
directly  to  DTRA  in  July  2014  and  December  2015. 

In  addition,  CRAFT  Tech  has  actively  participated  in  the  December  2015  Simulant  Chemistry 
Workshop  organized  by  DTRA.  This  was  an  excellent  opportunity  to  further  strengthen  the  link 
between  DTRA's  Basic  Research  portfolio  and  the  applications  side  of  DTRA's  activities. 
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